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Preface 


This  work  attempts  to  provide  the  solution  to  a problem 
first  identified  by  Major  J.  W.  Carl  of  the  Air  Force  In- 
stitute of  Technology,  School  of  Engineering.  The  basic 
problem  was  to  design  and  implement  the  software  necessary 
to  allow  two-dimensional  linear  filter  and  two-dimensional 
Discrete  Fourier  Transform  (DFT)  to  be  performed  on  an 
existing  minicomputer-based  image  processing  system.  With 
the  addition  of  the  routines  described  in  this  thesis, 
true  linear  filtering  of  images  is  now  available  on  the 

DICOMED  Image  Processing  System.  It  is  hoped  that  with  its 
filtering  capability,  and  with  its  64x64  DFT  capability 
that  this  new  image  analysis  software  will  indeed  provide 
a useful  analytical  tool. 

I would  like  to  take  this  brief  opportunity  to  thank 
the  people  without  whom  this  project  would  have  been  vir- 
tually impossible.  Thanks  to  Captain  Fred  Barney  for  his 
making  the  DICOMED  equipment  available  and  for  his  patient 
explanations  on  the  intracacier  of  same. 

A special  thanks  to  my  thesis  advisor.  Major  Carl  for 
his  saint-like  patience  and  profound  guidance  during  some 
very  trying  times.  And  lastly,  but  by  no  means  least, 
thanks  to  my  wife  Marla  for  her  typing  the  thesis  rough 
draft  and  for  her  inspirational  moral  support. 


Dennis  A.  McGaugh 
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Abstract 


A FORTRAN  program  is  presented  that  will  interactively 
prepare  digitized  images  (pictures)  for  experimental  analysi 
The  digital  preparation  includes  the  ability  to  Discrete 
Fourier  Transform  (DFT) , Inverse  Fourier  Transform  (IDFT) , 
threshold  image  power  spectral  density,  and  to  perform  two- 
dimensional  image  filtering.  Further,  a brief  discussion 
of  image  distortion  measurements,  optimal  image  size  for 
experimental  presentation,  and  image  aliasing  is  provided. 

A sample  of  typical  processed  image  results  is  also  supplied 
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AN  INTERACTIVE  ANALYSIS  PROGRAM  FOR 
A DIGITAL  IMAGE  PROCESSING  LABORATORY 


I.  Introduction 
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There  currently  exists,  located  in  the  Air  Force  Avionics 
Laboratory,  image  processing  hardware  and  software  capable 
of  digitizing  and  performing  limited  statistical  calculations 
of  pictorial  imagery.  There  also  exists  a requirement  to 
perform  certain  image  processing  with  this  equipment  in 
preparation  for  a visual  experiment.  However,  the  software 
needed  to  conduct  the  processing  did  not  formerly  exist. 

The  types  of  computations  needed  to  perform  this  processing 
and  analysis  are  well-known.  But  these  computations  needed 
to  be  developed  in  software  and  implemented  on  the  particular 
small  computer  system  that  is  available. 

The  image  processing  functions  that  are  needed  include: 
two-dimensional  linear  filtering,  calculation  of  a two-dimen- 
sional image  power  spectrum,  ability  to  specify  a threshold 
for  this  pov^er  spectrum,  and  the  Mean-Square-Error  (MSE) 
and  Bandwidth  (BW)  of  an  image  constrained  to  occupy  only 
those  frequencies  at  which  the  power  spectrum  is  above  a 
specified  threshold.  Additionally,  it  is  desired  that 
these  functions  be  interactively  controlled  and  specified 
to  the  program  during  run  time. 

The  problem  is  to  develop  and  provide  a software  pro- 
gram capable  of  accomplishing  the  above  functions,  and  still 
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be  compatible  V7ith  the  existing  equipment.  Additionally, 
part  of  this  thesis  work  is  to  determine  what  is  the  minimum 
size  of  picture  digitization  (e.g.  64x64,  128x128,  256x256) 
that  will  produce  pictures  suitable  for  the  experiment;  that 
is,  a picture  that  creates  an  image  of  no  less  quality  than 
30  lines  per  1°  of  visual  angle  subtended  during  viewing. 

The  picture  digitization  will  also  impact  on  the  computer 
resources  needed  to  perform  the  filtering.  The  experiment 
currently  calls  for  the  use  of  up  to  24  different  pictures. 
These  pictures  are  available  in  positive  or  negative  form, 
and  have  been  digitized  (2048x2048)  and  stored  on  magnetic 
tape.  The  sampling  resolution  may  also  cause  image  aliasing 
problems,  and  this  effect  must  also  be  discussed. 

The  material  presented  in  this  thesis  follows  roughly 
the  same  order  as  discussed  above.  Chapter  II  discusses  the 
constraining  environment,  its  effect  on  program  development, 
and  the  digital  domain  for  implementation  of  the  required 
functions.  Chapter  III  briefly  describes  the  laboratory 
hardware/software  and  the  final  project  software.  Chapter 

IV  presents  a sample  of  processed  image  results,  and  Chapter 

V closes  with  summary  remarks  and  mention  of  areas  for 
future  work.  The  appendix  contains  a program  user's  guide 
and  complete  program  documentation. 


II. 


Discussion 


During  the  development  of  the  image  analysis  program, 
several  important  questions  had  to  be  answered.  The  answers 
to  some  are  rather  straight-forward;  the  ansv/ers  to  others 
were  influenced  by  answers  to  previous  questions.  The 
answer  to  each  question,  and  its  impact  on  the  design  of 
the  image  analysis  program  is  presented  in  the  section  that 
follows'.  After  reading  these  sections,  it  is  hoped  that 
there  will  be  clear  understanding  of  the  program  constraint 
environment. 


f 
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A Question  of  Image  Resolution 

The  requirement  exists  to  conducb  visual  experiments 
with  pictorial  presentation  to  human  subjects.  The  images 
should  fall  within  the  subjects'  regions  of  foveal  vision; 
that  is,  the  images  should  not  exceed  a visual  angle  of  2 
degrees  during  viewing.  In  discussions  with  (then)  Captain 
Larry  Goble  of  the  US  Air  Force  Flight  Dynamics  Laboratory 
at  VJright-Patterson  AFB,  Ohio,  the  point  was  made  that  the 
industry  standard  for  images  on  graphics  and  television 
displays  is  a twelve-inch  square  picture  viewed  at  a distance 
of  thirty  inches  with  a resolution  of  no  less  than  sixty 
lines  per  inch.  It  was  also  decided  that  this  should  be  the 
resolution  standard  for  the  analysis  program.  After  some 
basic  geometrical  calculations,  this  works  out  to  be  30 
lines  of  resolution  per  1 visual  degree  (vis) . Additionally, 
the  experimental  displays  require  that  each  visual  presentation 
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consist  of  tv;o  different  pictures  placed  side-by-side.  Con- 
sidering the  overall  limit  of  2°  vis  vertically  and  hori- 
zontally, this  means  that  for  square  pictures,  each  picture 
can  only  subtend  1°  vis. 

The  ansv/er  to  the  question  of  image  resolution  is  that 
each  picture,  in  order  to  subtend  1®  vis,  must  be  1/2  inch 
X 1/2  inch  at  30  inches  distance  from  the  viewer,  and  each 
picture  must  have  a minimum  of  30  lines  of  resolution  both 
horizontally  and  vertically.  Imagery  of  this  quality  is 
obtainable  through  the  use  of  the  laboratory  DICOMED  image 
processing  equipment.  It  should  be  noted,  however,  that 
the  requirement  for  number  of  lines  of  resolution  increases 
when  considering  the  question  of  aliasing. 

A Question  cf  Aliasing  in  Digitized  Images 

The  DICOMED  image  processing  equipment  digitizes  a 
given  picture  v/ith  a sequential  horizontal  scanning  action. 

It  begins  at  the  top  of  the  picture  to  be  digitized  with 
a horizontal  scan  the  width  of  the  picture.  It  then  proceeds 
incrementally  down  the  picture  with  successive,  picture 
width,  horizontal  scans  until  it  stops  at  the  picture  bottom. 
Both  horizontal  and  vertical  directions  are  discretely 
sampled  with  an  adjustable  number  of  samples  per  picture 
widths  (independently  chosen  in  each  direction).  The  result 
of  this  scanning  process  is  the  mapping  of  the  entire  picture 
into  a MxN  matrix  of  M*N  discrete  pixel  values,  where  M and 
N are  both  powers  of  tv/o.  In  this  application,  M was  equal 


to  N.  For  the  DICOMED  machine,  the  value  of  each  pixel  can 
be  any  one  of  256  gray  scale  levels,  from  white  (377g)  to 
black  (OOOg) . When  referring  to  the  resolution  of  this  re- 
sulting digitized  picture,  it  is  said  to  have  N lines  of 
resolution,  with  a sampling  rate  of  N (where  N=2^) . 

The  answer  to  the  question  of  aliasing  can  be  found 
by  referring  to  the  Nyquist  sampling  theorm  (12:68).  Re- 
stated, it  simply  says:  to  eliirdnate  spectrum  aliasing  in 
the  discrete  domain,  sample  at  a rate  greater  than  or  equal 
to  twice  the  spectral  width.  Therefore,  for  the  previous 
case  of  an  MxN  sampled  image,  in  order  to  prevent  aliasing 
after  discrete  Fourier  transformations,  there  must  not  be 
any  spatial-frequency  content  above  E cycles  per  picture 
width.  The  experimental  requirement  is  for  a minimum  of 
30  lines  resolution  per  picture.  This  means  that  a digitized 
picture  with  this  resolution  can  display  up  to  30  full  os- 
cillations (cycles)  betv/eon  white  and  black  per  picture  in 
either  the  horizontal  or  vortical  direction.  Therefore, 

the  display  resolution  criterion  is  satisfied  if  > 30. 

2 

Since  a DFT/FFT  computer  algor ithim  (discussed  ]ater)  will 
be  used  to  transform  into  the  Fourier  domain,  and  powers 
of  tv/o  run  fastest  in  that  algorithim,  and  since  the  DICOMED 
image  processing  system  samples  at  powers  of  two,  a sampling 
rate  of  64  was  selected  as  the  minimum  convenient  value  that 
N can  have  and  still  meet  the  display  resolution  criteria. 
This  requires  a 64x64  digitization  of  each  imago  to  be  pro- 
cessed, and  insures  that,  at  this  digitization  level,  the 
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resolution  limits  of  human  observers  will  be  met  under  the 
cited  viewing  conditions. 

Now,  to  insure  no  image  aliasing,  it  must  be  guaranteed 
that  there  is  no  spatial  frequency,  in  the  image  being  trans 
formed,  above  32  cycles  per  picture  width.  To  provide  this 
guarantee,  all  64x64  digital  images  used  in  this  work  were 
derived  from  pictures  sampled  at  a rate  much  higher  than  64 
samples  per  line.  The  original  pictures  were  first  digi- 
tized into  2048x2043  images.  Then  by  considering  64x64 
blocks  of  32x32  pixels,  each  32x32  pixel  block  was  averaged 
into  one  pixel  value.  The  resultant  averaged  pixel  values 
produce  the  required  64x64  digital  image.  In  essence,  this 
procedure  performs  a two-dimensional,  digital  low-pass  fil- 
tering operation  on  a picture  that  may  contain  spatial  fre- 
quencies of  up  to  1024  cycles.  Although  this  type  of  low- 
pass  filtering  results  in  a form  for  the  filter  in  the 

X 

frequency  domain,  and  such  filters  are  not  particularly  good 
in  the  stop-band  of  the  filter  (considerable  side-lobe  struc 
ture)  , these  effects  v;ere  ignored  since  relatively  low 
spatial  frequencies  are  anticipated  in  the  original  pictures 
This  procedure  is  essentially  a digital  simulation  of  the 
low-pass  filtering  operation  necessary  to  avoid  aliasing. 


The  Question  of  Linear  Filtering 
The  functional  requirement 
spatial  convolution  of  the  form 
y(nj_,n2)  = x{n2,n2)  » h(nj^,n2) 


Implementa  tion 
for  linear  filtering 


(10M40) 
^ J 


mp=-()C 


mn=-C>0 


implies 


h(nj^-mj^,  n2-m2) 
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where  y(nj^,n2)  is  the  output  image  result  of  the  x(nj^,n2)  J 

input  image  that  has  been  filtered  by  a linear  filter  with 

I 

an  impulse  response  of  h(n^,n2).  As  noted  by  Rabiner  and  ^ 

Gold  (Ref  10) , this  calculation  in  the  spatial  domain  with 

images  is  not  usually  soluable  in  closed  form  except  for  the 

most  trivial  filters.  However,  a more  tractable  form  of 

this  convolution  equation  occurs  in  the  frequency  domain. 

It  is  well-known  (Ref  8)  that  the  convolution  equation  takes 
% 

the  form: 

. y(ki,k2)  = X(ki,k2)  * HCk^,  k2) 
where  Y(kj^,k2),  X(kj^,k2),  and  II(k2^,k2)  are  the  Fourier  trans- 
form coefficients  of  the  image  output,  image  input,  and 
image  filter  impulse  response  respectively.  Now,  the  filter 
computations  become  readily  manageable  and  almost  trivial., 

They  simply  involve  multiplying  each  Fourier  transform 
coefficient  of  the  input  image  with  each  respective  Fourier 
transform  coefficient  of  the  impulse  response  of  the  filter. 

The  result  is  then  inverse  Fourier  transformed  back  into  the 
spacial  domain  as  the  resultant  linearly  filtered  output 
image . 

The  transfer  function  furnished  for  the  linear 
filtering  process  is  given  as  a one  dimensional  equation. 

Before  it  can  be  used  to  filter  an  image,  however,  it  must  be 
converted  into  two-dimensional  form.  It  has  been  shown  by 
Huang  (Ref  6)  that  a good  circularly  symmetric  two-dimensional 
linear  filter  can  be  obtained  from  a good  one-dimensional 
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linear  filter  via  the  relation; 

H(k3^,k2)  = 

A 

where  H is  the  one-dimensional  filter  impulse  response  at 
the  appropriate  frequency  values  of  kq  and  k2.  For  the 
situation  involving  the  previously  cited  64x64  image  matrices, 
this  requires  Fourier  domain  frequency  values  (k3^,k2),  along 
both  the  real  and  imaginary  axis,  ranging  from  -32  cycles 
to  +32  cycles.  Due  to  the  circular  symmetry  of  the  filter, 
only  one  32x32  quadrant  need  be  calculated  from  the  equation. 
Once  calculated,  it  is  appropriately  copied  into  the  other 
three  quadrants  of  the  filter  matrix.  A more  complete  ex- 
planation of  this  operation  appears  in  Appendix  A. 

The  Question  of  Operating  in  Tv/o  Dimensions 

The  most  obvious  problem  of  processing  in  two-dimensions 

is  simple:  it  takes  longer,  and  uses  more  computer  memory. 

For  the  one-dimensional  case,  the  discrete  Fourier  transform 

2 

of  a sequence  of  N samples  takes  N complex  addition  and 
multiplications;  in  two-dimensions,  it  takes  operations 
by  brute  force.  For  the  one-dimensional  case,  the  filtering 
process  in  the  Fourier  domain  requires  2«N  complex  coeffic- 
ients to  be  stored  in  computer  memory;  in  two-dimens  ons 
2 

2«N  complex  coefficients  must  be  stored.  To  help  reduce 
these  increased  requirements  associated  with  two-dimensional 
processing,  certain  processing  modifications  were  made. 

First  was  the  implementation  of  the  well-known  (Ref  8) 
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Fast  Fourier  Transform  (FFT)  technique.  This  reduces  the 
number  of  addition  and  multiplications  from  (64^=16,77,216) 
to  2N*Nlog2N  (128*64-6  = 49,152),  a reduction  of  over  10 
to  1.  The  program  implementation  of  the  FFT  algorithim  is 
via  the  FORTRAN  subroutine  called  FOURT.  It  is  a program 
written  by  Norman  Brenner  from  the  basic  program  by  Charles 
Rader  of  MIT  Lincoln  Laboratory.  It  is  a very  versatile 
program,  and  offers  more  options  than  are  required  for  the 
particular  image  processing  application.  But  its  very 
versatility  was  the  reason  it  was  selected,  as  it  will  allow 
calculation  of  odd  sizes  of  FFT's  for  images  other  than  the 
present  64x64  size.  A more  complete  discussion  of  FOURT 
appears  with  the  program  listing  in  the  appendix. 

The  second  modification  was  to  use  disk  files  to  store 
the  images  and  complex  FFT  coefficients.  In  the  existing 
IPS/DICOMED  image  processing  systems  each  complex  number 
occupies  4 words  of  storage.  For  one  64x64  FFT  coefficient 
matrix  to  be  entirely  loaded  into  computer  memory  it  would 
require  16K  words.  If,  as  in  a filtering  operation,  two 
matrices  were  to  be  used  at  once,  then  32K  words  would  be 
needed;  the  existing  minicomputer  core  storage  is  only  28K 
words.  A technique  suggested  by  Oppenhein  (Ref  8)  to  reduce 
the  requirement  for  two-dimensional  core  storage  is  to  per- 
form the  two-dimensional  FFT  as  a series  of  one-dimensional 
FFT's.  Further,  by  first  storing  the  images  on  random  access 
disk  files  then  reading  in  only  one  row  or  column  from  each 
appropriate  disk  file,  the  FFT  operates  on  only  those  two 
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pieces  of  data  at  one  time,  and  the  core  required  to  operate 
the  program  is  dropped  to  a realizeable  size.  Although  the 
processing  scheme  involves  huge  amounts  of  computer  I/O, 
the  time  to  compute  and  store  a 64x64  point  complex  FFT  is 
less  than  4 minutes,  and  the  time  to  accomplish  filtering 
(matrix  multiplication)  is  less  than  2 minutes.  Addition- 
ally, there  is  an  option  to  have  the  FFT  coefficients  nor- 
malized by  dividing  all  coefficients  by  the  zero-frequency 
or  DC  term,  or  to  have  the  FFT  coefficients  remain  with 
their  unnormalized  values. 

^ Question  of  Power  Spectrum  Thresholding 

Once  the  digitized  image  is  in  the  Fourier  domain,  the 
calculation  of  the  power  spectrum  is  straight-forward.  The 
power  spectral  density  of  an  image  is  related  to  the  squares 
of  the  absolute  value  of  the  Fourier  coefficients  and  is  given 
by: 

^(u,v)  = E ||x(u,v)| 

where  the  expectation  is  over  an  ensemble  of  image  transforms. 
The  program  implementation  for  estimating  the  power  spectrum 
of  a set  of  images  is  obvious  from  this  equation,  and  a 
discussion  can  be  found  in  Appendix  A,  Operating  Procedures, 
5c.  The  analysis  j-J^ogram  provides  the  capability  to  average 
up  to  24  different  image  spectrums.  These  24  images  may  be 
the  same  ones  that  are  mentioned  in  the  introduction  to  this 
paper.  It  is  on  this  resultant  averaged  power  spectrum  that 


the  thresholding  function  is  expected  to  be  used. 

The  thresholding  function  was  originally  requested  to 
be  used  with  only  the  power  spectrum.  However,  to  provide 
more  flexibility,  it  is  implemented  in  a manner  to  allow 
its  use  with  any  specified  64x64  matrix  of  complex  data 
(i.e.  FFT  coefficients,  or  power  coefficients).  By  using 
the  thresholding  function,  the  analysis  program  user  is 
able  to  interactively  specify  a 64x64  complex  valued  matrix, 
and  to  enter  a real  valued  number  against  which  the  absolute 
value  of  every  complex  point  is  compared.  The  result  of 
this  thresholding  operation  is  the  creation  of  a 64x64  matrix 
which  contains  (in  complex  format)  the  value  one  for  every 
point  above  the  threshold  and  a complex  zero  for  every  point 
below  the  threshold.  Both  the  original  matrix  and  the  thresh- 
old matrix  are  then  available  for  further  processing. 

Since  the  thresholding  is  done  in  the  Fourier  domain, 
values  that  are  linearly  related  (Ref  1:138-141)  to  spatial 
bandwidth  (BW)  and  mean-square-error  (MSB)  are  easily  cal- 
culated. The  EW  value  is  simply  the  number  of  points  above 
threshold  in  the  Fourier  domain.  In  the  spatial  domain, 
these  above-threshold  values  correspond  to  the  image  pixels 
that  would  be  transmitted  over  a digital  channel  with  relative 
bandv/idth  BW.  The  MSB  value  is  simply  the  sum  of  the  power 
values  below  threshold.  In  the  spatial  domain,  the  MSB 
corresponds  to  the  error  of  the  image  transmitted  over  a 
digital  channel,  with  bandwidth  BW  described  above. 
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The  Question  of  Images  and  Image  File  Generation 

Originally,  the  only  requirement  for  image  output  was 
to  generate  a completely  processed  image  ready  for  experi-  ^ 

mental  testing.  However,  in  an  effort  to  make  the  analysis 
program  applicable  for  more  general  use,  the  ability  to  dis- 
play intermediate  results  is  also  provided.  A 64x64  image 
file  displaying  the  results  for  each  of  the  following 
functional  operations  is  prepared  and  stored  in  DICOMED/IPS 
format  on  disk;  (1)  the  reduction  of  a 2048x2048  image  to 
a 64x64  image;  (2)  the  two-dimensional  FFT  calculation;  (3) 
the  power  spectrum  calculation;  and  (4)  a fully  processed 
image  (either  thresholded,  linearly  filtered,  or  FFT'd) . 

Since  the  images  are  in  DICOMED/IPS  format,  post-pro- 
cessing  with  the  IPS  system  is  trivially  implemented  (Ref 
7) . A sample  image  file  format  appears  in  the  appendix. 
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III.  Laboratory  Hardware  and  Software 

The  imagery  equipment  used  for  this  work  is  located  in 

I 

the  Avionics  Laboratory  on  VJright-Patterson  APB,  Ohio.  A 
detailed  description  of  the  equipment  is  included  in  Refer- 
ence 7.  Briefly,  it  consists  of  a PDP-11/45  minicomputer 
with  28K  of  core  memory  and  many  peripherals.  Included  are: 
one  RP^3  diskpack  disk  drive,  tv;o  RK^5  diskcar tridge  disk 
drives,  one  MTUljS  tape  drive,  one  132  column  line  printer, 
a Conrac  video  monitor,  and  one  Tektronix  4000  series 
graphics  terminal.  Additionally,  the  system  is  interfaced 
to  a DICOMED  image  scanner /recorder . With  this  equipment, 
the  capability  exists  to  perform  up  to  a 2048x2048  point 
digitization  of  pictures.  As  the  DICOMED  digitizes  the 
picture,  it  is  stored  on  magnetic  tape  for  follow-on  pro- 
cessing by  the  PDP-11/45  Image  Processing  System  (IPS)  . A 
description  of  the  IPS  software  is  r' . o included  in  Refer- 
ence 7 . ► 

A full  description  of  the  software  written  for  this 
thesis  appears  in  the  appendix.  It  is  written  for  inter-  j, 

active  use  and  provides  the  following  capabilities: 

1.  Image  reduction  from  2048x2048  to  64x64. 

2.  FFT  of  a 64x64  image. 

3.  Calculation  of  .]^i(u,v)  for  one  image,  or  an  average 

A 

|)(u,v)  of  a standard  set  of  24  images. 

4.  Ability  to  threshold  a 64x64  image  as  previously 
discussed. 

5.  Ability  to  generate  a 64x64  filter  from  a user 
supplied  impulse  response  equation. 
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6.  Inverse  FFT  of  Fourier  coefficients. 


) 

I 

l 

f 

By  running  this  program,  any  one  of  the  above  processes  can  ^ 

produce  an  image  file  that  is  stored  on  the  RP^13  disk  for  , 

display  by  the  IPS  system,  or  hardcopy  picture  generation  by 
the  DICOMED  system. 


/ 
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IV.  Results 


Figure  1 presents  some  samples  of  output  from  the 
digital  image  preparation  program.  It  should  be  noted 
that  Figure  1 v/as  generated  with  an  arbitrary  high-pass 
filter  whoso  transfer  function  has  the  form: 


H{f)  = 0.5(1  + ln(f/10  + 0.1),  for  0 f ^ 22 


Also  note,  that  the  slight  graininess  of  the  64x64  pictures 
is  due  to  a blow-up  factor  of  two  (4  new  points  for  each 
original  point) . This  was  done  solely  to  provide  larger 
pictures  for  Figure  1. 


Figure  1.  Sample  Images 
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V.  Summary  and  Recommendations 


This  thesis  has  presented  a computer  program  to  digit- 
ally process  images,  and  in  this  case  the  processing  is 
motivated  by  a desire  to  generate  stimuli  for  a particular 
psychometric  experiment.  Additionally,  as  can  be  verified 
by  reading  the  program  description  in  the  appendix,  a gen- 
eral purpose  routine  is  now  made  available  to  tailor  digital 
processing  of  pictures  for  a large  number  of  image  filtering 
based  applications. 

Areas  for  future  work  would  include  modifying  the  pro- 
gram to  process  pictures  other  than  the  64x64  size,  and  an 
on-line  interactive  tie-in  between  the  IPS  system  and  the 
processing  program. 
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User's  Guide  for  RAZMA.TAZ 


Introduction 

This  program  allows  a user  to  perform  the  following 
functions  on  IPS/DICOMED  images: 

1.  A reduction  of  a 2048x2048  image  file  to  a 64x64 
image  file. 

2.  A tv;o-dimensional  fast-fourier  transform  (FFT)  of 

a 64x64  image  with  either  normalized  or  unnormalized 
coefficients . 

3.  Calculation  of  a two-dimensional  power  spectrum  of 
one  64x64  image  file,  or  the  two-dimensional  average 
power  of  24,  64x64  images. 

4.  Specification  of  a high-pass  threshold  on  either 
the  FFT  file  or  the  power  file. 

5.  Generation  of  a two-dimensional  filter  (FFT  domain) 
from  a user  specified  impulse  response  equation. 

6.  Multiplication  of  any  threshold  or  filter  files 
with  any  image  or  FFT  or  power  files. 

7.  Calculation  of  the  two-dimensional  inverse  FFT 
transform  on  any  modified  FFT  file. 

Reference  Material 

1 . The  DOS/DATCH  Handbook:  Monitor  Operating  System 

Version  09 . Digital  Equipment  Corporation,  Maynard, 


'■  'fK'. 
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DEC-ll-ODBllA-A-D,  April  1974  . 


Mass . 

2 . linage  Processing  SoftVviarc  Conversion , Vo  1 . I & II. 
Development  Center,  Air  Force  Systems  Command, 
Grifiss  AFB,  New  York  13441.  February  1977. 
(ADA038136  & ADA038137) . 

Operating  Procedures 

This  gives  the  ordinary  routine  for  running  RAZMA.TAZ. 
Specialized  uses  and  procedures  are  covered  under  sub-heading 
Miscellaneous  Operations. 

1 . Program  materials  you  will  need  include: 

User's  Guide. 

IPS  diskeartridge  (RK05)  (AAT  Lab-blue  label) . 

"Scratch"  diskeartiridge  (RK^5)  (AAT  Lab-#6) . 

ROLLIN  magtape  with  RAZMA.TAZ  program  and  24  picture 
data  sets. 

IPS/DJCOMED  PDP-1445  computer  and  associated  hard- 
ware (AAT  lab) . 

2.  How  to  'power-up': 

a.  Place  PDP-11/45  console  HALT  switch  down,  turn  key 
of  console  clock-wise ^ to  the  'on'  position. 

b.  Load  IPS  diskpack  in  RK^5  unit  #0.  Load  'scratch' 
diskpack  in  RKj35  unit  #1.  Do  not  use  'write  protect' 
switch.  VJait  until  'on-cylinder'  light  comes  on. 

c.  Mount  ROLLIN  tape  on  tape  drive.  This  can  be 
omitted  if  program  and  data  files  arc  already  on 
tlie  scratch  diskeartridge  or  the  RP0  disk  drive. 

d.  Turn  on  RAMTEK  graphics  display  box. 

e.  Slowly  turn  on  line  printer,  then  toggle  'master 
clear',  and  toggle  'on-line'  switches.  This  can 
bo  omitted  if  no  printer  output  is  expected  (nor- 
mally none  is) . 
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f.  Open  door  to  RPj3  disk  drive  to  see  if  diskpack  #1 
(DOS)  is  installed.  If  not,  retrieve  it  from  the 
lower  cabinet  and  install  it.  Turn  on  power  switch 
and  wait  until  yellow  light  comes  on. 

g.  Turn  on  Tektronix  4014  terminal  and  place  it  in 
'local'  mode.  Type  'ESC  key,  then  ':'  key.  Type 
in  that  order,  and  only  once  each.  This  reduces 
the  display  front  size.  Place  terminal  'on-line'. 

h.  Turn  on  CONRAC  video  display. 

i.  Return  to  the  PDP-11/45  console  SWR.  Set  in  000776g 
('one'  is  up)  and  toggle  'load  ADR'  switch.  Set 

in  000777g  and  toggle  'DEP'  switch.  Move  switch 
from  'HALT'  to  'ENABLE'  and  hit  the  'RUN'  sv.’itch. 

j.  Place  the  thumbwheel  of  the  AED8000  to  the  'MAINTEN- 
ANCE' display  and  to  ADDR  CODE  = 776732.  Push  the 
red  IPL  switch  once.  The  LED's  above  the  display 
should  begin  flashing.  If  not,  go  back  and  try 
again . 

k.  Return  to  the  Tektronix  terminal  and  type  RK/i  Do 
not  type  a carriage  return.  Enter  the  date  and  tim.e 
in  the  format  shown. 

l.  Login  with  $ LO  11,11. 

m.  Type  $ RUN  DKl;  RAZMA.TAZ,  or  $ RUN  DP0;  RAZMA.TAZ. 

n.  Follow  the  yellow  brick  road. 

3.  How  to  'power-down': 

a.  To  end  RAZMA.TAZ,  type  the  'CTRL'  key  and  the  'C' 
key  simultaneously,  followed  by  'K'  and  'I'.  This 
exits  to  the  DOS  monitor.  This  may  also  be  used 
to  stop  the  program  at  any  time.  Next,  type  'FI' 
to  log-out  from  the  IPS/DOS  system. 

b.  Place  the  terminal  in  'LOCAL'  mode,  give  the  screen 
display  several  RESET' s,  and  turn  off  the  terminal. 

c.  Turn  off  the  CONRAC  video  monitor. 

d.  Return  to  the  PDP-11/45  console  and  hit  the  'HALT' 
switch. 

e.  Toggle  the  switches  on  both  RK^5  diskcartridge  drives 
to  tlie  UNLOAD  position.  The  doors  will  remain 
locked  until  the  'LOAD'  light  comes  on.  This  will 
take  about  a minute.  Then  the  cartridges  can  be 
removed . 
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f.  Meanwhile,  the  following  steps  may  be  performed. 


g.  Toggle  the  Line  printer  to  'OFF  LINE'  and  turn  it 
off. 

h.  Turn  off  the  RAMTEK  graphics  display  box. 

i.  Push  the  red  IPL  switch  on  the  AED8000  once.  Push, 
then  release,  the  yellow  power  switch  on  the  RP^f 
disk  drive.  The  yellow  light  will  slowly  go  out. 
After  it  goes  out  the  diskpack  can  be  removed  if 
desired . 

j.  Unload/dismount  the  ROLLIN  magtape,  if  it  was  loaded 
on  the  magtape  drive,  and  turn  off  the  magtape  drive. 

k.  Return  to  the  RKj2I5  diskcartridge  drives  and  remove 
the  diskcartridges . 

l.  Turn  the  key  on  the  PDP-11-45  console  fully  counter- 
clockwise to  the  off  position. 

4.  Running  the  RAZMA.TAZ.  The  actual  operation  of  the 
program  is  straight  forward.  A menu  of  functions  is 
presented  first  and  the  user  selects  the  desired 
function.  All  prompts  are  self-explanatory.  The 
format  for  all  ASSIGN  commands  is  given  prior  to 
each  request  for  user  assignnient . Any  error  in 
syntax  will  be  caught  by  the  monitor  system.  However, 
be  very  careful  to  note  the  device  unit  number  used 
in  the  command  format.  Failure  to  use  the  same 
number  when  typing  in  the  command  can  destroy  certain 
necessary  stored  program  files  along  with  the  user's 
specified  input  file.  Should  the  wrong  number  be 
typed,  simply  retype  the  command  before  typing  'CO'. 

If  'CO'  has  already  been  typed,  do  a 'CNTRL/C,  'KI' 
immediately.  This  closes  all  files  and  returns 
control  to  the  DOS  monitor.  By  running  FILDMP  or 


PIP  (see  DOS/BATCK  handbook),  the  user's  files  can 
be  verified  and/or  restored.  For  verification  of 
program  files,  see  Special  Procedures.  After  file 
verification  or  restoration,  RAZMA.TAZ  can  be 
started  normally. 

Other  than  prompts  and  asking  for  file  assing- 
ments  there  should  be  no  other  messages,  except  when 
using  the  FFT,  POK’ER,  and  IFFT  functions.  In  these 
cases,  an  error  message  will  be  printed  about  an 
overflow  on  a real  to  integer  conversion.  This  is 
normal.  It  occurs  during  the  generation  of  each 
function's  image  file,  and  refers  to  the  DC  term  of 
the  FFT  coefficients.  It  does  not  affect  any  in 
ternal  processing  calculations.  The  error  merely 
means  that  the  value  of  the  DC  term  is  greater  than 
16,384  and  exceeds  the  allowable  256  level  gray 
scale  of  the  DICOMED  image  processor.  The  DC  pixel 
value  will  be  zero  (black  pixel) , in  the  upper  left 
hand  corner  of  any  picture  created  from  the  image 
file . 

Errors  of  any  other  type  should  not  occur.  If 
they  do,  it  is  most  likely  due  to  processing  a bad 
data  from  an  input  file.  As  before,  use  FILDMPT  or 
PIP  to  check  input  files  for  correct  data. 

5.  RAZMA.TAZ  functions.  This  briefly  describes  what 
each  function  does  and  what  files  are  used  or  gener- 
ated . 
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a.  REDUCTION.  Takes  a user  specified  2048x2048  image 
file  and  averages  it's  pixels  into  a 64x64  image 
file.  The  routine  takes  a 32x32  pixel  block  from 
the  specified  file  and  averages  the  entire  1024 
points  into  one  point  of  the  64x64  image.  The  out- 
put file  is  called  AVG.IMG. 

b.  FFT.  Performs  a two-dimensional  fast  fourier  trans- 
form (FFT)  on  either  a user  specified  image  file  or 
the  AVG.IMG.  file.  The  resultant  coefficients  can 
be  either  normalized  (FFTCOF.NRM  file)  or  unnormal- 
ized (FFTCOF.UNN  file)  as  requested  by  the  user.  The 
routine  performs  the  FFT  by  using  the  subroutine 
FOURT,  first  on  each  image  row,  and  then  on  each 
image  column.  An  image  file  of  the  resultant  FFT 
coefficients  is  also  made  (FFT. IMG).  The  DC  term 

is  in  the  upper  left-hand  corner.  This  image  file 
may  be  accessed  in  the  normal  manner  under  the  IPS/ 
DICOMED  system. 

c.  POWER.  Uses  the  FFT  coefficients  to  generate  a 
two-dim.ensional  power  spectrum.  The  input  file 
can  be  a single  user  specified  FFT  file,  or  FFTCOF. 
NRM,  or  FFTCOF.UNN.  Additionally,  the  function 
will  calculate  the  average  power  spectrum  of  the 

24  storage  picture  coefficients.  The  power  cal- 
culation is  simply  the  squaring  of  the  FFT  coeffi- 
cients. The  resultant  output  file  is  PWRFFT.COF. 

An  image  file  is  also  created  (PWRFTT.IMG)  for  use 
on  the  IPS/DICOMED  system. 

d.  THRESHOLD.  This  function  looks  a user  specified 
coefficient  file,  or  PWRFFT.COF,  or  FFTCOF.NRM,  or 
FFTCOF. UMM  and  creates  a file  that  when  used  with 
the  MULT  function  will  pass  only  coefficients  above 
a certain  threshold  value.  This  value  is  entered, 
interactively,  by  the  user.  In  addition,  the 
routine  counts  the  number  of  points  above  the 
threshold  and  calculates  a "bandwidth"  number;  i.e., 
PTS. ABOVE  f 4096  = BW.  Also,  if  the  coefficients 
were  POWER  coefficients  and  normalized,  it  finds 
the  mean-square-error  of  the  thresholded  power 
spectrum;  i.e.,  the  amount  of  power  below  the 
threshold  setting.  The  output  file  contains  only 

a mask  of  I's  and  J3 ' s at  the  appropriate  location 
and  is  called  THRHLD.FIL. 

e.  FILTER.  This  function  takes  a user  specified  filter 
transfer  function  equation  and  computes  a two-dimen- 
sional mask  file  which  can  be  used  in  the  MULT 
program.  The  value  of  the  mask  points  are  based  on 
the  radial  distance  from  the  upper  left-hand  corner 
of  a 32x32  square  (Fig  2) . 


Figure  2.  User  Specified  Filter  Impulse  Response 


The  equation  is  specified  on  line  515  of  the  RAZMA.TAZ 
program  listing.  To  change  the  equation,  the  pro- 
gram must  be  edited  and  reassembled  (see  Miscellan- 
eous Operations,  b) . The  resultant  output  file  is 
FILTER.COF. 

f.  MULT.  This  functional  routine  takes  one  user  spec- 
ified coefficient  file,  or  THRHOLD.COF,  or  FILTER. 

COF  and  multiplies  it  point-by-point  times  another 
user  specified  coefficient  file,  or  FFTCOF.NRM,  or 

FFTCOF.UNN,  or  PWRFFT.COF.  It's  output  file  is  ^ 

MULT.DAT. 

g.  IFFT.  This  function  performs  an  inverse  FFT  on 
the  MULT.DAT  file.  The  user  cannot  directly  specify 
input  to  this  file.  But  he  can  however,  input  in- 
directly. By  giving  the  THRESHOLD  function  a 
threshold  of  zero,  a mask  of  all  one's  will  be 
generated.  Then,  using  this  mask  with  MULT  and 
specifying  the  other  MULT  input  file,  the  input  to 
IFFT  is  assured.  The  output  is  an  image  file 
named  PROCSD.IMG.  As  with  all  other  image  files, 
this  is  operable  on  the  IPS/DICOMEu  system. 

6.  FLOWCHART.  The  flowchart  in  Figure  3 shows  the 

processing  flow  and  options  available  in  RAZMA.TAZ. 

f 
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Special  Procedures 

1.  NO  FILE!  After  the  $LO  and  $RUN  the  monitor  ansv;ers 
with  NO  FILE!  When  this  occurs,  run  the  DOS/DATCII 
PIP  program  (see  the  handbook)  and  see  if  RAZMA.TAZ 

is  on  the  disk  you  specified  in  the  $RUN  command.  ’ 

If  not,  check  the  other  disk.  Remember,  you  must 

be  logged-in  under  UIC=11,11.  If  not  on  either 

disk,  use  the  ROLLIN  tape  and  ROLLIN  program  (in 

handbook)  to  restore  the  RK^5  diskcartridge  on 

unit  #1  (DKl:).  Then  start  POWER-UP  procedures 

2.  UNABLE  TO  OPEN  FILE.  This  can  occur  any  time 
during  the  running  of  the  RAZMA.TAZ  program.  This 
usually  means  that  an  expected  input  file  (spec- 
ified by  ASSIGN)  or  one  of  the  24  standard  picture 
files  does  either  not  exist  on  DP^:  or  two  files 
exist  with  the  same  name.  In  either  case,  run 
PIP  to  get  a listing  of  all  files  on  DP)3:  to  see 
if  all  required  files  are  there,  and  perform  file 
manifulation  as  appropriate.  If  file  cannot  be 
found  on  DP)3:  or  DKl:,  it  may  be  necessary  to 
ROLLIN  the  magtape  as  in  the  above  step  and  start 
again.  A minimum  of  the  following  files  are  re- 
quired: 
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DISCCTOIW  DltJi  e : 

11. n 3 

lS-S£P-7a 

ARTSl  .cor 

64C 

13-£EP-7fl 

<233> 

AftTya  .COF 

64C 

13-fEP-78 

<e33> 

ARTyi  .COF 

G4C 

13-CEP-78 

<83  3> 

ASTy^  .COf 

e4C 

13-f.E.P-78 

<e33> 

APCl  .COF 

64C 

13-SEP-7a 

<83  3> 

APC2  .COF 

64C 

13-£EP-7fi 

<:e33> 

APC3  .COF 

64C 

13-5EP-7ft 

<233; 

Af>f4  .COF 

e4C 

l3-feP-78 

<233> 

CORTICI.COF 

64C 

n-6EP-V8 

''2337 

cwetKa.coF 

(>4C 

J3-SEP-78 

<2337 

CWTK3.C0F 

64C 

J3-SEP-7S 

<233- 

CVRTK^.COF 

640 

13-SEP-78 

<233/ 

OPNTKl.CCif 

64C 

i3-scp-7a 

(233/ 

0PNTIC2.C0F 

64C 

13-6EP-7B 

<2337 

OPNTKa.COF 

64C 

13-SEP-78 

'233.- 

OPNTK4.COF 

64C 

13-SCP-7a 

<2337 

fHERl  .COF 

64C 

13-SEP-78 

(233: 

£KLP2  .COF 

64C 

l3-?EP-7a 

(233; 

SHtPJ  .COf 

64C 

13- SEP -78 

(2337 

SHtP4  .COF 

640 

13-5EP-78 

'833. 

nCOAl  .COF 

64C 

13-SEP-78 

<2337 

BCOAC  .COF 

64C 

13-5EP-78 

■:233; 

P5CA3  .COf 

64C 

13-fEP-78 

.2  3 3; 

n6«A4  .COF 

64C 

13-SEP-78 

<233: 

ARTyi  .me 

lec 

13-SEP-76 
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APTya  , 

mo 
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<2335 

ARTVa  , 

me 
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<2337 

ARTy4  , 
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13-SEP-78 
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IMG 
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IMG 
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PASHA  , 

FTM 
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IMN 

TOIL  FlLCSl 

M 

i 

I 
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3.  EXPONENT  OVERFLOW  or  OVERFLOW  DURING  REAL  TO  INTEGER 
CONVERSION.  This  can  occur  during  the  operation 

of  all  functions  except  FILTERING  or  THRESHOLDING. 

It  is  the  result  of  processing  with  invalid  data  in 
the  assigned  input  files,  or  bad  data  from  one  of 
the  24  stored  image  files.  In  either  case,  the 
contents  of  files  can  be  inspected  using  the  DOS/ 
BATCH  program  FILDMP  (see  handbook) . The  organi- 
zation and  contents  of  typical  image  and  coefficient 
files  are  shown  in  Appendix  C.  If  bad  data  is 
suspected  in  any  of  the  24  noted  image  files,  the 
best  procedure  is  to  use  the  DOS/BATCH  program  PIP 
to  restore  the  files;  an  example  of  how  to  do  this 
is  shovm  in  the  next  section. 

4.  UNABLE  TO  READ  SPECIFIED  FILE  FORMAT  or  END  OF  FILE- 
END  OF  MEDIUM.  This  occurs  when  the  RAZMA.TAZ  pro- 
gram encounters  files  with  im.proper  file  organization 
or  length.  Again,  as  above,  use  FILDMP  or  PIP  to 
verify  and  transfer  files. 


Miscellaneous  Operations 

1.  Generating  Files.  There  are  three  categories  of 
files  associated  with  the  RAZMA.TAZ  program. 

a.  Programming  operating  files.  These  are  files 
generated  by  the  program  for  it's  own  use.  A dis- 
cussion of  these  files  is  found  in  the  comments  of 
the  program  listing. 

b.  Twenty-four  standard  picture  file  sets.  Each  pic- 
ture set  consists  of  an  image  file  and  a normalized 
FFT  coefficient  file.  These  picture  sets  have 
already  been  generated  and  are  stored  for  use  by 
the  program  in  the  POWER  function.  A copy  of  these 
files  is  on  the  ROLLIN  magtape. 

c.  User  impact  files.  These  are  files  specified  by  the 
user  with  the  ASSIGN  command.  The  file  will  be 
either  the  image  or  the  FFT  coefficient  type. 

The  origin  of  the  image  files  can  be  from  the  IPS/ 
DICOMED  system  (2048x2048  or  64x64)  or  the  RAZMA.TAZ 
program  (64x64).  Directions  for  creating  image 
files  under  the  IPS/DICOMED  system  can  be  found  in 
Volume  1 of  Reference  2.  Creating  a 64x64  image 
file  under  RAZMA.TAZ  can  be  accomplished  by  using 
the  reduction  function,  and  then  exiting  to  use 
PIP  to  rename  the  newly  generated  AVG.IMG  file  to 
whatever  name  the  user  desires.  To  save  these 
files  on  magtape,  PIP  and  ROLLIN  are  used.  For 
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each  file,  the  PIP  command  is:  #DK1:  name.extn/ 
CO<DP0:  name.extn.  Then,  ROLLIN  is  used  (DOS/ 

BATCH  handbook)  to  copy  DKl:  and  create  a new 
ROLLIN  tape. 

The  origin  of  the  FFT  coefficient  files  must  be 
from  running  RAZMA.TAZ.  The  file  is  simply  created 
by  running  the  FFT  function,  and  then  exiting  to 
use  PIP  to  rename  the  newly  created  FFTCOF.NRM  or 
FFTCOF.UNN  file  to  whatever  name  the  user  desires. 
The  file  can  then  be  saved  in  the  same  manner  as 
the  image  files. 

2.  Restoring  files.  For  RAZMA.TAZ  to  operate,  all 
files  must  be  on  DP0: . If  they  are  not  (a  minimum 
listing  appears  in  Special  Procedures),  or  for 
some  reasons  the  files  become  invalid,  the  files 
must  be  restored.  A permanent  copy  of  all  files 
is  on  the  ROLLIN  magtape.  This  tape  can  be  loaded 
onto  DKl:  (DOS/BATCH  handbook),  and  then  the 
appropriate  file  transferred  under  PIP  to  DP0: . 

The  proper  PIP  command  is:  #DP0:  name . extn/CO < 

DKl;  name.extn. 

3.  Changing  RAZMA.TAZ.  Best  advice:  DON'T.  However, 
for  the  brave  and/or  foolhardy,  a bubble  chart  and 
program  listing  is  given  in  Appendix  B.  Once  the 
changes  have  been  made  under  EDIT  (see  handbook) , 
the  proper  Fortran  assembly  must  be  made.  The 
proper  command  after  $RUN  FORTRAN  is;  #DKl;  RAZMA, 
LP:<DK1;  edit  file  name,  extn  / ON;  or,  #DP0; 

RAZMA,  LP;<DK1;  edit  file  name.  extn/ON.  (Assuming 
EDIT  was  done  on  DKl;).  This  is  important  because 
RAZMA.TAZ  uses  single  word  integer  variables.  When 
the  assembly  is  done,  the  linker  command  after 

$RUN  LINK  is:  #DK1;  RAZMA.TAZ < DKl^  RAZMA/CC,  DK0: 
FTNLIB/L/E;  or;  #DP0:  RAZMA . TAZ  < DKl ; RAZMA/CC, 

DK)3:  FTNLIB/L/E  as  appropriate. 
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User ' s Guide 
Program  Bubble  Chart 
Program  Listing 


NAME  IMG 


IMAGE 

AVRG 

(100) 


Figure  4.  Program  Bubble  Chart 
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PROGRAM  RAZMA.TAZ 

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 


THIS  PROGRAM  has  THE  FOLLOWING  DISK  (DP0|)  RANDOM  ACCESS  FILE 
ASSIGNMENTS! 

AVG  ,1MG  • OUTPUT  OF  REDUCE#  INPUT  TO  FPT 
SCRATC.DAT  • USED  IN  FFT 

FFTCOF.NRM  « OUTPUT  OF  FFT#  INPUT  TO  POWER 
(FFTCOF.UNN)  input  to  THRESHOLD 

INPUT  TO  MULT 

FFT  ,IHG  - OUTPUT  OF  FFT 

,,COF  • INPUT  TO  POWER 

PWRFFT.COF  • OUTPUT  OF  POWER#  INPUT  TO  THRESHOLD 
PWRFFT.IMG  - OUTPUT  OF  POWER 

THRHLO.FIL  • OUTPUT  OF  THRESHOLD#  INPUT  TO  MULT 
FILTER, COF  • OUTPUT  OF  FILTER#  INPUT  TO  MULT 
MULT  ,DAT  • OUTPUT  OF  MULT#  USED  IN  IFFT 
PR0C8D,IMC  • OUTPUT  OF  IFFT 

THE  FOLLOWING  DEVICES  ARE  INTERACTIVELY  ASSIGNED  TO  THESE 
FILES  DURING  THE  RUNNING  OF  THE  INDICATED  MODULES! 


FILE 

1 

2 

3 

4 

5 6 

7 

BIG  ,IMG 

REDUCE 

AVG  ,IMG 

FFTCOF,*— 

MULT 

THRE8 

POWER 

FFT 

powfft,cof 

FILTER, COF 

MULT 

THRE8 

THRMLD,COF 

MULT 

ANY  I/O  REFERENCE  MADE  TO  DEVICE  NUMBERS  AFTER  THE  DEVICES  HAVE 
BEEN  INTERACTIVELY  ASSIGNED  DURING  THE  MODULES  SPECIFIED  ABOVE, 
WILL  OPERATE  ON  THE  FILE  DATA  EQUATED  WITH  THAT  DEVICE  NUMBER, 


CUHPLEX  0ATA(6R)#0ATBUF(6R}#DCVAL#N0FILS#M8E 
BYTE  HEADER (64# 8) # BUFFER (2045) #BBUF (2) #BCOL( 126)# NULLS (6a# 6) 
integer  INPTR#AVGPTR#DAT3PT#0AT2PT#FFTPTR»DATTPT#DAT8PT 
INTEGER  8UM#BUF#0FF8ET#C0L(64)#UN1T#C0FPTR#C0FIM6#AV6FFT 
integer  THR8PT#0UTBW#BW, IFRR,DAT1PT#DAT5PT#DAT4PT# A8GD7#PWR0N 
EQUIVALENCE  (BCOL ( 1) # COL ( 1 ) ) # (BBUF ( 1 ) # BUF) 

DATA  HEAOCR/0#0#79#56#64#0#64#50S*e/ 

DATA  NULL8/512*0/ 

1 WRITE  (6#5) 

WRITE  (6#6} 

WRITE  (6#U) 

2 WRITE  (6#7) 

WRITE  (6#e) 
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I 

< 

) 


WRITE  (6»9) 

WRITE  (6»10) 

. WRITE  (6,11) 

5 PORWAT(*  GREETINGS  AND  SOLUTIONS  *) 

6 FORMATC'  SPEAK  TO  ME  ABOUT  IMAGE  INPUT  *) 

7 FORHATf'  TTPEl  D FOR  BIGGIE  /20R8X204«/  ') 

8 FORMAT!*  J FOR  JUNIOR  /6aX6«/  *) 

9 FORMAT!'  S FOR  STANDARD  SET  OF  29  ') 

10  FORMAT!'  N FOR  NONE  ') 

11  FORMAT!'  ') 

read  !6,12)  INBJSN 

12  F0RMAT!A) 

IF!1NBJ8N,E0,1HJ)  GOTO  200 
IF!INBJSN,EQ(1HB}  GOTO  100 
iFdNBJSN.EQ.lHN)  GOTO  30 
IF!INBJSN,NE,1H8)  GOTO  2 

20  WRITE  !6,25) 

WRITE  !6,26) 

WRITE  !6,27) 

22  WRITE  !6,28) 

WRITE  !6,11) 

25  FORMAT!'  THE  STANDARD  SET  OF  29  PICTURES  IS  ALREADY  STORED  ') 

26  FORMAT!'  IN  IMAGE  FILES  fc  IN  NORMALIZED  FFT  COEFF  FILES  ON  DP0|  ') 

27  FORMAT!'  IF  WISH  TO  PROCEED  TO  DO  POWER,  FILTER,  OR  THRESHOLD,  ') 

28  FORMAT!'  TYPE  C TO  CONTINUE,  TYPE  S TO  STOP,  ') 

READ  (6,12)  INCS 

IF!INC8,EQ,1HS}  GOTO  900 
IF!INC8,NE,1HC)  GOTO  22 
30  WRITE  !6,35) 

WRITE  !6,36) 

WRITE  (6,37) 

WRITE  (6,38) 

WRITE  (6,11) 


35 

FORMAT!' 

TYPEl  1 

TO 

DO 

FFT  ') 

36 

FORMAT!' 

2 

TO 

DO 

POWER  ') 

37 

FORMAT!' 

i 

TO 

DO 

THRESHOLDING  ') 

38 

FORMAT!' 

9 

TO 

DO 

FILTERING  ') 

READ  (6,12)  INHERE 
ir!IWHERE,EQ,lHl)  GOTO  200 
IF!IWHERE,EQ,1H2)  GOTO  300 
IF!IWHERE,EQ,1H3}  GOTO  900 
IF!IWHERE,NE,1H9)  GOTO  30 
GOTO  500 

C tk*****  REDUCTION  MODULE  ***** 

100  WRITE  (6,105) 

WRITE  (6,107) 

WRITE  (6,108) 

WRITE  (6,109) 

105  FORMAT('  THE  INPUT  IMAGE  FILE  MUST  BE  SPECIFIED,  ') 

107  FORMAT!'  IF  NOT  YET  A88ICNE0»  PLEASE  DO  SO  NOW,  THE  FORM  I8|  ') 

108  FORMAT!'  SA8  DP0|NAME,EXT, 7 ') 

109  FORMAT!'  THEN  TYPE  CO  TO  CONTINUE,  ') 

A86D7»'Y' 

PAUSE 

DEFINE  FILE  7(9097, 512, U, INPTR) 
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I 

1 

I 


CALL  SETFIL(8f 'AVC,IMGMERR,'OP*,0) 
define  file  8(80, 32,U, AVGPTR) 

NRITE  (6,110) 

110  FORHAT(*  AVERAGE  WORKING  *) 

INPTRa2 
AVCPTRal 
DO  115  J«l,8 

115  WRITE  (S'AVCPTR)  (HE ADER ( I , J) , I«1 , 64) 

DO  116  Jal>6 

116  WRITE  (8'AVCPTR)  ( NULLS ( I , J) , I* 1 » 64) 

DO  170  1*1,64 
DO  120  K«l,64 
120  COL(K)*0 

DO  150  J*l,32 

READ  (7MNPTR)  (8UFFER(N)  , N*1 , 1024) 

READ  (7'INPTR)  ( BUFFER ( N) , N* 1 025, 2048) 

DO  140  K*l,64 

0FFSET*(K-l)wi2 

SUH«0 

DO  130  L«lf32 
INOXaLtOFFSET 
B8UE( l)*BUFFER(INOX) 

SUMaSUHf BUF 

130  CONTINUE 

C0L(K)aSUH/326C0L(K) 

140  CONTINUE 

150  continue 

DO  160  Kal,64 
C0L(K)*C0L(K)/32 
160  CONTINUE 

WRITE  (8'AVGPTR)  (BCOL ( N) , Nal , 128, 2) 

170  CONTINUE 

ENDFILE  7 
FNDFZLE  8 

190  WRITE  (6^195) 

WRITE  (6,11) 

192  WRITE  (6,196) 

WRITE  (6,11) 

195  FORMAT!*  DONE  IMAGE  AVERAGING  *) 

196  FORMAT!*  TYPE  C TO  CONTINUE,  TYPE  8 TO  STOP  *) 

READ  (6,12)  INCS 

IF(1NC8.CQ,1H8)  GOTO  900 
iPdNCS.NE.lHC)  GOTO  192 
GOTO  30 

C mmti*  FFT  MODULE  ***** 

200  call  8ETF1L(8,*AVC,IMC*,IERR,*DP*,0) 

DEFINE  FILE  8(80, 32, U, AVGPTR) 

IF(A8GD7.CQ.1HY)  GOTO  215 
WRITE  (6,205) 

WRITE  (6»207) 

WRITE  (6,208) 

WRITE  (6*209) 

205  FORMAT!*  THE  INPUT  IMAGE  FILE  MUST  BE  SPECIFIED,  *) 

207  FORMAT!*  IF  NOT  YET  ASSIGNED*  PLEASE  DO  SO  NOW,  THE  FORM  ISI  *) 
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2H8  FORMAT('  SAS  DP0 I MAME , EXT , 7 ') 

' 289  FORMAT('  THEN  TYPE  CO  TO  CONTINUE,  *) 

. PAUSE 

1 DEFINE  FILE  7(80,32,U,AVGPTR) 

> 210  PRITE  (6,215) 

. WRITE  (6,216) 

WHITE  (6,11) 

I 215  FORMAT('  TYPEl  N FOR  NORMALIZED  FFT  COEFF  •) 

' 216  FORMAT(»  U FOR  UNNORMAUZED  COEFF  *) 

I READ  (6,12)  NORMUN 

DEFINE  file  3(64, 256, U,DAT3PT) 

* CALL  8ETFIL(3,*8CRATC,DAT*,IERR,'DP»,0) 

. DEFINE  file  2(64, 256, U, DAT2PT) 

I IF(NORMUN,tO,lHN)  CALL  SETF1L(2, 'FF TCOF  , NRMMERR,  »OP ' , 0) 

i IF(NORMUN,EO, IHU)  CALL  SE7F IL (2, 'FFTCOF , UNN» , lERR , ’DP * , 0) 

. 1F(N0RHUN,NE, IHU, AND, NORMUN, NE.IHN)  GOTO  210 

217  DEFINE  FILE  1 (80, 32, U, FF TPTR) 

CALL  StTFIL(  1, 'FFT,1mg»,  lERR, 'DP',«') 

' WRITE  (6,219) 

I 219  FORMAT(*  FFT  WORKING  ') 

AVCPTR«17 

DAT2PT«1 

DAT3PT*! 

FFTPTR»1 

. SCALE»l,0/64,0 

DO  220  J«1,0 

220  WRITE  (I'FFTPTR)  (HEADER ( I , J) , I«l , 64) 

I DO  225  J«t,8 

' 225  WRITE  (1*FFTPTR)  ( NULLS ( I , J) , I« 1 , 64) 

I DO  230  Kb1,64 

230  COL(K)s0 

! DO  250  lsl,6a 

; IF(ASG07,EQ,1HY)  REAO(e*AVGPTR)  (BCOL ( N) , Nsl , 128, 2) 

I IF(AS0D7,NE,1HY)  REAO(7*AV6PTR)  (BCOL ( N) , N« 1 , 128, 2) 

i DO  240  Jal,64 

> DATA(J)sCOL(J} 

’ 240  CONTINUE 

) CALL  FOURT(DATA, 64,1, *1,1,0) 

> WRITE  (3'0AT3PT)  (DAT A ( N) , N«l , 64) 

I 250  CONTINUE 

DO  290  K«l,64 

f DATlPTai 

* DO  260  I«l,64 

\ READ  (3»0AT3PT)  (DATBUF ( N) , N«1 , 64) 

> OATA(I)ROATBUF(K) 

> 260  continue 

' CALL  F0URT(DATA,64,1,. 1,1,0) 

I IF(NORMUN,EO,IHU)  GOTO  270 

* 8CALEa256, 0*64,0 

) J«1 

IF(K*1)  265,265,266 
! 265  DCVALaDATA(l) 

i J«J*1 

I 266  DO  267  IsJ,64 

» OATA(I)aDATA(I)/OCVAL 


I 

t 


— — . jr  . 


I 

) 

1 
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continue 

write  (2'DAT2PT)  (DATA(N),Nsl,6a) 

DO  260  I«l,6<l 

COl,(I)«fCABS(DATA(n))wSCALE 

continue 

WRITE  (I'FFTPTR)  ( BCOt.  ( N) , Na  1 , 1 28, 2) 

CONTINUE 
ENDFILE  I 
ENOFIUE  2 
ENOFILE  3 
ENDFILE  7 
ENOFILE  8 
WRITE  (6,296) 

WRITE  (6,11) 

WRITE  (6,297) 

WRITE  (6,11) 

FORMAT!*  DONE  FFT  *) 

FORMAT!*  TYPE  C TO  CONTINUE,  TYPE  8 TO  STOP  ') 

READ  (6,12)  INCS 
IF(INC8,Ea,lHS)  GOTO  900 
IF(INCS,Nt,lHC)  GOTO  295 
WRITE  (6,299) 

WRITE  (6,11) 

FORMAT!*  TYPE  P FOR  POWER,  F FOR  FILTER,  T FOR  THRESHOLD  *) 

READ  (6,12)  IPFT 
IFdPFT.EO,  IHF)  GOTO  500 
1F(IPFT,EQ,IMT)  GOTO  «00 
IF(IPFT,NE,IHP)  GOTO  298 
GOTO  300 


C 

300 


i 

• 301 

I 

302 

303 

I 

i 30a 


I 


I 


I 

i 305 

. 306 

307 

I 

I 306 


*****  POWER  MODULE  ***** 

CALL  8ETFIL(8, *PWRFFT , I MG* , lERR, »DP*,0) 

DEFINE  FILE  8(80, 32,U,DAT8PT) 

DATSPTai 
DO  301  Jsl,8 

WRITE (e*DAT8PT)  ( HEADER ( I , J) , lai , 64) 

DO  302  J»l,e 

WRITEC8*DAr8PT)  ( NULLS ( I , J) , lal , 64) 

WRITE(6,304) 

WRITE  (6,11) 

FORMAT!*  TYPEl  0,  FOR  ONLY  ONf  FFT  INPUT)  S,  FOR  STD  24  FFT  COfcFF  *) 

READ  (6,12)  INOS 

PWRONa*Y* 

IF(IN08,E0,tHS)  GOTO  312 
IF(IN08,NEtlH0)  GOTO  303 

IF(NORMUN,EO,lHN)  CALL  8ETFIL(5, *FFTCOF , NRM* , lERR, 'DP* , 0) 
IF(N0RMUN,E0,1HU)  CALL  SETF IL (5, *FFTCOF , UNN* , lERR, *DP * , 0) 
IF(NORMUN,EQ, 1HN«0R«N0RMUN,CQ,1HU)  GOTO  308 
WRITE  (6,305) 

WRITE  (6,306) 

WRITE  (6,307) 

FORMAT!*  FFT  INPUT  MUST  BE  SPECIFIED,  PLEASE  DO  SO  NOW,  *) 

FORMAT!'  THE  FORM  I8l  SA8  OP0|NAME,EXT# 5 *) 

FORMAT!*  THEN  TYPE  CO  TO  CONTINUE,  •) 

PAUSE 

define  FILE  5(64, 256, U,DAT5PT) 
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3P9 


310 


311 


312 


313 


314 


315 

316 

317 

318 

319 

320 

321 


322 

323 


i 324 
325 


I 


CALL  SCTFILl  W*PWRFFT,COFMERR,»DP',0) 
DFFINE  FILE  1 (64, 256, U, DAT IPT) 
kxRITE  (6,309) 

FORHAT('  SQUARING  FFT/S/  ') 

DAT1PT«1 

DAT5PT*! 

DO  311  I«lf64 

READ  (5»DAT5PT)  ( DAT  A ( N) , Na 1 , 64 ) 

DO  310  Jsl,64 

DATA(  J)«DATA(  J)<.DATA(  J) 

COL(J)*(CABS(OATA(J))) 

CONTINUE 

WRITE  (I'DATIPT)  ( DAT A ( N) , N= I , 64) 

WRITE(8»0AT8PT)  ( BCOL ( N) , N» 1 , 128, 2) 

CONTINUE 

ENDFILE  1 

ENOFILE  5 

ENDFILE  8 

GOTO  390 

WRITE  (6,309) 

CALL  8ETFIL(5,'PWRFFT,C0F%URR, 

DEFINE  FILE  5 ( 64 , 256, U, D AT5PT ) 

SCALEb64,0**3 

N0F1LS«(24, 0,24,0) 

DO  513  I*l»64 
DATBUF(I)*(0, 0,0,0) 

CONTINUE 
DAT5PT«1 
DO  314  I8l,64 

WRITE  (5'DAT5PT)  (OATBUF ( N) , N«1 , 64) 

CONTINUE 

UNIT»5 

DO  360  1*1,24 

DAT5PT»1 

COFPTR*! 

IF(UNIT*4)  319,315,317 
DO  316  J*l,4 
ENDFILE  J 
DO  318  J*l,4 

DEFINE  file  J ( 64, 256, U, COFPTR)  . 
UNIT*0 
UNIT«UNlTfl 
IF(I«1)  322,32U322 

CALL  SETFILdr 'APC1,COF',URR,'DP',0) 
CALL  8ETFIL(2,'APC2,COFMERR,'DP»,0) 
CALL  SETFIL(3»'APC3,COF*,IERR,»DP»,0) 
CALL  8ETFIL(4,*APC4,C0F', 1ERR,'DP',0) 
IF(I-5)  335,323f324 

CALL  8ETFIL(1,'ARTY1,C0FMERR,  »DP*»0) 
CALL  8ETFIL(2,»ARTY2,COFMERR,»OP*,0) 
CALL  StTFILC3»'ARTY3,C0FMERR,  •OP*,0) 
CALL  SETFIL(R»'ARTY4,COF',IERR,*DP»,0) 
IF(I-9)  335,325^326 

CALL  SETFILd,  'CVRTKl,COF',IERR,  'OP»,0) 
CALL  8ETFlL(2r 'CVRTK2,C0FMERR,  'ORdB) 


1 


; 
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CAUL  SETFILC3#*CVRTK3,C0F', lERR, 'DP',0) 

CALL  SETFILC^r  'CVRTKiJ.COF  S lERR,  'DP*,0) 

526  IF(1-13)  335,527,328  _ 

327  CALL  SETFIL(l,'OPNTKl,COF'»IERP»'DP'»0) 

CALL  StTFIL(2,'0PNTK2,C0F',IERR, »DP',0) 

CALL  8E7FIL(3, '0PNTK3,C0F*, lERR, 'DP',®) 

CALL  SETFILCtl/'OPNTKfl.COFSlERR,  'DP',0) 

328  IF(I*17)  335,329,330 

329  CALL  SETFILd, 'M60AI,COF',  lERR,  »DP',05 
CALL  SETFIL(2, 'H60A2,COr', lEHR, »OP',0) 

CALL  SETFIL(3,'M60A3,COFMERR,»DP',0) 

CALL  SETFIL(9,  •M60A/|,COF',  lERR, 'DPS0) 

330  lF(l*?n  335,331,335 

331  CALL  SETFIL(l,'SHF.Rl,COFMERR,»DP',0) 

CALL  SETFILC2, 'SHER2,C0F', lERP, 'DP',0) 

CALL  SETFIL(3,'SHER3,COF»,IERR,'DP',0) 

CALL  SETFIU(«,'SHER9,C0F', lERR, 'DP»,0) 

335  CONTINUE 

00  350  J«l,6(i 

REA0(5*0AT5PT)  ( DATBUF ( N) , N* I , 6«) 

READ(UNIT'COFPTR)  ( 0 AT A ( N) , N« 1 , 6«) 

DO  337  K«l,6« 

537  OATA(K)«DATA(K)*DATA(K) 

DO  390  Kcl,69 

DATBUF(K)*DATBUF(K)+DATA(K) 

' 590  CONTINUE 

DAT5PTaDAT5PT«l 

WRITE(5*DAT5PT)  ( D ATBUF ( N) , Na I , 69) 

350  CONTINUE 

360  CONTINUE 

365  0AT5PT*1 

. DO  380  Js1,69 

READC5*DAT5PT)  ( DAT  A ( N) , N« i , 69 ) 

> DO  370  J»l,69 

' 0ATA(J)a0ATA( J)/NOFILS 

I COL(J)«(CAB8(OATA(J)))*SCALE 

370  CONTINUE 

DAT5PT»DAT5PT*l 

WRITE(5*DAT5PT)  (DAT A ( N) , N* I , 69) 

WRITE (e»DAT8PT)  ( BCOL ( N) , Na 1 , 1 28 , 2) 

. 380  CONTINUE 

DO  385  lBl,5 
385  ENOFILE  I 

ENDFILE  8 

' 390  WRITE  (6,395) 

• WRITE  (6,11) 

395  FORMAT('  DONE  POWER  ') 

; GOTO  900 

c *****  threshold  module  ***** 

I 900  WRITE  (6,905) 

I WRITE(6,11) 

; 905  FORMAT(*  TYPE  N TO  MARE  NEW  THRESHOLD  FILE, TYPE  0 TO  USE  OLD  FILE 

• READ  (6,12)  INNO 
IF(INN0,E0, IHO)  GOTO  600 

i iFdNNO.NE,  IHN)  GOTO  900 
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410 

WRITE  (6,415) 

WRITE  (6,416) 

WRITF  (6,11) 

415 

FORMAT!'  TYPEl 

V FOR 

416 

FORMAT!' 

A FOR 

VANILLA  PLAIN  FFT  COEFF  ') 

averaged/souared  fft  input  n 


READ  (6»12)  INVA 

IKINVA.EQ.IHA)  CALL  SETF  IL  ( 4, ‘PWREFT , COP  ' , I ERR , ' OP  ' , 0) 

IF ( I NV A, EQ, IHV.AND.NORMUN, EQ, IHN) 

1C  ALL  StTP  IL(4,  »PFTCOF,NRM%  lERR,  'DP',0) 

IF(INVA,EQ, IHV, ANO.NORMUN.EQ, IHU) 

I CALL  SETFIL(4,  'FFTCOF  , UNNS  ILRR » 'DP  ' , 0 ) 

IFdNVA, NE.IHA, AND, INVA. NE.IHV)  GOTO  <410 

lF(NORMUN,tQ, IHN.OR.NORMUN.EQ, IHU.OR.PwRON.EQ, IHY)  GOTO  4i0 
WRITE  (6»434) 

WRITE  (6#425) 

WRITE  (6»426) 

WRITE  (6,427) 

WRITE  (6,426) 

424  FORMAT('  THE  FILE  TO  BE  THRESHOLDED  MUST  BE  8PECIPIED,  ') 

425  FORMAT('  IF  NOT  YET  ASSIGNED,  PLEASE  DO  SO  NOW,  ') 

426  FORMAT!'  THE  FORM  IS|  ') 

427  FORMAT!'  SAS  DP0 1 NAME , EX T , 4 ') 

428  FORMAT!'  THEN  TYPE  CO  TO  CONTINUE  ') 

PAUSE 

430  define  file  4 ( 64 , 256, U, DAT4PT) 

call  SETFILd,  'THRHLD.FIL',  TERR,  'DP',0) 

DEFINE  FILE  1 ( 64 , 256, U, THRSPT ) 

WRITE  (6,431) 

432  WRITE  (6,433) 

WRITE  (6,11) 

431  FORMAT!'  SPECIFY  A REAL  VALUE  FOR  THE  THRESHOLD  OF  ') 

433  FORMAT!'  THE  FORM  * 0,0000  /UP  TO  15  TRAILING  DIGITS/  ') 

READ  (6,435,ERR«432)  THRES 

435  F0RMAT(Fi8, 15) 

WRITE(6,439) 

WRITE  (6,11) 

439  FORMAT!'  THRESHOLD  WORKING  ') 


DAT4PT«1 

THRSPT»1 

BW30 

MSE»(0, 0,0,0) 

DO  470  1*1,64 

READ  (4'DAT4PT)  ( D AT  A ( N) , N» 1 , 64 ) 

DO  460  J8l,64 
DATBUF! J)*( 1,0,0.0) 
BELOW*CABS(OATA(J)  )-THRF.S 

450  IF(BELOW)  451,453,453 

451  MSE*MSE+DATA( J) 

0UTBW*0UTBW+1 
DATBUF(J)«(0. 0,0.0) 

453  CONTINUE 

460  CONTINUE 

WRITE d' THRSPT)  ( D AT BUF ( N) , N*1 , 64) 
470  continue 
EN0F2LE  1 
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ENDFILt  a 

BM*4096*OUTbW 

ABSMSt«CAHS(MSt) 

FRACBWBFUOAT(8W)/i|096,0 

NRITE(6,a80)  IMRES 
WRITE  (6,48?)  ABSMSE 
WRITE  (6,484)  BW,FRACBW 
WRITE  (6,11) 

480  FORMAT!'  THIS  THRHLO.FIL  USED  A REAL  THRESHOLD  OF  ',F10,15) 

482  FORMAT!'  THE  ABSOLUTE  MSE  WAS  s ',F18,1B) 

484  FORMAT('  AND  THE  BANDWIDTH  wAS*  ',l4,'/4096  OR  ',F12,R) 

490  WRITE  (6,495) 

WRITE  (6,11) 

491  WRITE  (6,496) 

WRITE  (6,11) 

49b  FORMAT!'  END  THRESHOLD  ') 

496  FORMAT!'  TYPE  C TO  CONT INUE, TYPE  S TO  STOP  ') 

READ  (6,12)  INCS 
IF(  INCS.EQ,  IHS)  GOTO  900 
IFdNCS.NE,  IHC)  GOTO  491 
GOTO  600 

C *****  FILTER  MODULE  ***** 

500  WRITE  (6,505) 

WRITE  (6,506) 

WRITE  (6,11) 

505  FORMAT!'  TYPEl  C TO  USE  CURRENT  FILTER  EfiUATION  FOR  FILE  *) 

506  FORMAT!'  0 TO  USE  AN  OLD  FILTER  FILE  ') 

READ  (6,12)  INCO 

IF (INCO.EQ, IHO)  GOTO  600 
IF (INCO, NE, IHC)  GOTO  500 
510  WRITE  (6,512) 

512  FORMAT!'  MAKING  FILTER  ') 

DEFINE  FILE  1 (64, 256, U, OAT IPT ) 

CALL  8ETFIL(l,'FILTER,C0F', lERR, 'DP',0) 

KAPA1B0 

KAPA2«0 

SIGMA1»0, 

SIGMA2a0, 

DAT1PT«1 
1F0LD*34 
DO  530  m*1,64 

MMbM 

DO  520  N»l,64 

NNaN 

IF (M,CE, IFOLD)  MMB66-M 

IbMM«1 

IF(N,CE, IFOLD)  NNB66-N 

JaNN»l 

8QlJal**2+J**? 

RADIALaSORT(SQIJ) 

515  DATA(N}b0,5*(  1.0* ALOG (RADI  AL*0|U1«0) ) 

520  CONTINUE 

WRITE  (I'DATIPT)  (DAT A ( N) , Na 1 , 64) 

530  CONTINUE 


r 


T 


j 

I 

1 
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ENOFiLt:  1 * 

WRITE  lb,bH5) 

WRITE  (6,U) 

540  WRITE  (6»5a6) 

WRITE  (6,11) 

545  PORHAT('  DONE  FILTER  ») 

546  FORMATt*  TYPE  C TO  CONTINUE,  TYPE  S TO  STOP  ') 

READ  (6,12)  INCS 

IFdNCS.EQ,  IMS)  GOTO  900 
IFdNCS.NE.lHC)  GOTO  540 
GOTO  600 


C *****  matrix  MULT  MODULE  ***** 

600  IFdNNO.EQ,  IMN)  CALL  SE  TF IL  d , 'THRHLD , F IL  ' , lERR,  » DP  • , 0 ) 
IF(INCO,Ea,lHC)  CALL  St TF IL (1 , ' F ILTER , COF • , lERR, • DP  * , 0) 
TF(1NnO,NE,1HO,AND,1NCO,NE,1HO)  GOTO  620 
IFdNNO.EQ, IHO)  WRITE  (6,605) 

IFdNCO.EQ,  IHO)  WRITE  (6,606) 

605  F0RMAT('  the  old  threshold  file  must  be  SPECIFIED,  d 

606  FORMAT(»  THE  OLD  FILTER  FILE  MUST  BE  SPECIFIED,  d 

WRITE  (6,607) 

WRITE  (6,606) 

WRITE  (6,609) 

WRITE  (6,610) 

607  FORMAT('  IF  ALREADY  ASSIGNED,  TYPE  CO  TO  CONTINUE,  d 

608  FORMAT!'  IF  NOT  YET  ASSIGNED,  PLEASE  DC  SO  NOW,  THE  FORM  ISl  d 

609  FORMAT!'  SAS  DP0 1 NAME , EXT , I ') 

610  FORMAT!'  THEN  TYPE  CO  TO  CONTINUE,  ') 

PAUSE 

620  DEFINE  FILE  I (64, 256, U, DAT IPT ) 

IF(NORMUN,tO,lHU)  CALL  SETF 1 L ( 3, 'FFTCOF , UNN' , lERR , ' DP ' , 0) 
IF(N0RMUN,EQ,1HN)  CALL  SETF I L C 3, 'FFTCOF , NRM' , lERR , ' DP ' , 0) 

WRITE  (6,625) 

WRITE  (6,626) 

WRITE  (6,627) 

WRITE  (6,628) 

WRITE  (6,629) 

625  FORMAT!'  IF  FFT  COEFF  WERE  JUST  NOW  MADE,  THE  FILE  IS  ') 

626  FORMAT!'  ALREADY  ASSIGNED,  TYPE  CO  TO  CONTINUE,  ') 

627  FORMAT!'  IF  NOT  YET  ASSIGNED,  PLEASE  DO  SO  NOW,  THE  FORM  IS|  d 

628  FORMAT!'  SAS  OP0 | NAME , EXT , 3 ') 

629  FORMAT!'  THEN  TYPE  CO  TO  CONTINUE,  ') 

PAUSE 

define  FILt  3(64, 256, U,DAT3PT) 

IF(NORMUN,EQ, 1HU,0R,N0RMUN,EQ,  IHN)  GOTO  640 

630  WRITE  (6,635) 

WRITE  (6,11) 

635  FORMAT!'  TYPE  N FOR  NORM  FFT  COEFF,  TYPE  U FOR  UNNORMALIZED  COEFF  ') 
READ  (6,12)  NORMUN 

IF(NORMUN,NE, 1HU,AND,N0RMUN,NE,1HN)  GOTO  630 
640  CALL  8ETFIL(2,'MUiT,DAT',IERR,'DP',0) 

DEFINE  FILE  2 (64, 256, U, DAT2PT) 

WRITE  (6,645) 

645  FORMAT!'  MULTIPLYING  ') 

DATIPT«1 


\ 


W 
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DAT2PT«1 

0AT3PT»l 

DO  680  Isl»6a 

READ  (5'OATJPT)  (DATA(N),n* 

READ  (I'OATIPT)  (DATBUKN), 

IF(NORHuN,tO, IHU)  GOTO  665 


1>64) 

N«l»6a) 


» 

K«1 

j 

650»650»655 

> 

650 

WRITE  (6^651) 

r 

651 

FORMAT!'  AH»HA,,,,THINK  YOU  CAN  FOOL 

i 

DCVAL«OATA(l) 

1 

KaK4l 

) 

655 

DO  660  JBK»64 

1 

DATA(J)aDATA( J)*DCVAL 

> 

660 

CONTINUE 

i 

665 

DO  670  Jat,64 

1 

DATA(J)aDATA( J)aDATBUF( J) 

) 

670 

CONTINUE 

) 

write  !?'0AT2PT)  (DAT A (N) , Na l , 64) 

r 

680 

CONTINUE 

) 

ENDFILE  1 

) 

ENDFILE  2 

) 

ENDFILE  3 

1 

WRITE  (6»695) 

> 

WRITE  (6,11) 

5 

695 

FORMAT!'  done  mult  ') 

» 

GOTO  700 

C 

*****  inverse  fft  module  ***** 

> 

700 

CALL  8ETFIL(2,'MULT,DAT', IERR,'OP',0) 

> 

define  file  2(64, 256, U,DAT2PT) 

r 

CALL  8ETFIL(8,'PR0C80,IMG',IERR,'DP', 

) 

define  file  8(80,32,U,DAT8PT) 

> 

WRITE  (6,730) 

5 

730 

FORMAT  ('  IFFT  WORKING  ') 

1 

0AT2PTal 

> 

DATBPTal 

i 

SCALEal, 0/4096.0 

» 

DO  740  Jai.8 

> 

WRITE  (8'0AT8PT)  (HEADER ( I , J) , lai , 64) 

> 

740 

CONTINUE 

t 

DO  745  Jal,e 

) 

WRITE  (8'DAT8PT)  ( NULLS ( I, J) , lal , 64) 

) 

745 

CONTINUE 

) 

DO  750  lal,64 

1 

READ  (2'DAT2PT)  (DAT A ( N) , Nal , 64) 

> 

CALL  FOURT(DATA,64,l,l,i,0) 

1 

DAT2PTaOAT2PT*l 

» 

WRITE  (2'DAT2PT)  (DAT A ( N) , Nal , 64) 

> 

750 

CONTINUE 

) 

DO  700  K8l,64 

f 

DAT2PTal 

) 

DO  760  181.64 

> 

READ  (2'0AT2PT)  (DATBUF (N) , Nal , 64) 

5 

DATA(I)aDATBUF(K) 

MOTHER  NATURE  EHl  ') 
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CONTINUE 

CAUL  F0URT(DATA,64, t, 1» l»0) 

DO  770  I»l»64 

COL(n«(WEAL<DATA(I)))*SCALt 
770  CONTINUE 

NRITE  (8'DAT8PT)  (BCOt(N) , Nsl, 128, 2) 
780  CONTINUE 

ENDFILE  2 
ENDFUE  8 
WRITE  (6,795) 

795  FORMAT!'  DONE  IFFT  ') 

900  END 


ROUTINES  CALLEDi 

8ETFIL,  FOURT  , CABS  , FLOAT  , SORT  , ALOG  , REAL 


OPTIONS  ■/0N,/0Pl3 

BLOCK  LENGTH 

MAIN,  9645  (045532)* 

**COMPILER  CORE** 

PHASE  USED  FREE 

DECLARATIVES  00963  14349 
EXECUTABLES  03343  11969 
ASSEMBLY  04731  15221 
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SUBROUTINE  FOURT ( D AT A, NN, NDIM, IS 16N, I FORM, WORK) 

THE  COOLfcVTUKEY  FAST  FOURIER  TRANSFORM  IN  USASI  BASIC  FORTRAN 


TRANSF0RM(K1,K2, ,,,)  « SUH( OAT A ( J i , J2, , , . ) *E XP ( IS 1CN*2»PI*SURT ( • 1) 
*CCJl*J)*(Kl»n/NN(t)  + (J2-l)*(Ka-l)/NN(2)^, ,,))),  SUMMED  FOR  ALL 
Jl,  Kl  BETWEEN  I AND  NNC1)»  J2,  K2  BETWEEN  1 AND  NN(2),  ETC, 

THERE  IS  NO  limit  TO  THE  NUMBER  OF  SUBSCRIPTS,  DATA  IS  A 
multidimensional  complex  array  whose  real  and  IMAGINARY 
PARTS  ARE  ADJACENT  IN  STORAGE,  SUCH  AS  FORTRAN  IV  PLACES  THEM, 

IF  ALL  imaginary  PARTS  ARE  ZERO  (DATA  ARE  DISGUISED  REAL),  SET 
IFORM  TO  ZERO  TO  CUT  THE  RUNNING  TIME  BY  UP  TO  FORTY  PERCENT, 
OTHERWISE,  IFORM  a ♦!,  THE  LENGTHS  OF  ALL  DIMENSIONS  ARE 
STORED  IN  ARRAY  NN,  OF  LENGTH  NDIM,  THEY  MAY  BE  ANY  POSITIVE 
INTEGERS,  THO  THE  PROGRAM  RUNS  PASTER  ON  COMPOSITE  INTEGERS,  AMO 
ESPECIALLY  FAST  ON  NUMBERS  RICH  IN  FACTORS  OF  TWO,  ISICN  IS  ♦! 

OR  -I,  IF  A •!  TRANSFORM  IS  FOLLOWED  BY  A ONE  <OR  A 
BY  A .1)  THE  ORIGINAL  DATA  REAPPEAR,  MULTIPLIED  BY  NTOT  (aNN(l)w 
NN(2)a,,,),  TRANSFORM  VALUES  ARE  ALWAYS  COMPLEX,  AND  ARE  RETURNED 
IN  ARRAY  DATA,  REPLACING  THE  INPUT,  IN  ADDITION,  IF  ALL 
DIMENSIONS  ARE  NOT  POWERS  OF  TWO,  ARRAY  WORK  MUST  BE  SUPPLIED, 
complex  of  length  EUUAL  to  the  largest  non  2a«rK  DIMENSION, 
OTHERWISE,  REPLACE  WORK  BY  ZERO  IN  THE  CALLING  SEQUENCE, 

NORMAL  FORTRAN  DATA  ORDERING  IS  EXPECTED,  FIRST  SUBSCRIPT  VARYING 
fastest,  ALL  SUBSCRIPTS  BEGIN  AT  ONE. 

RUNNING  TIME  IS  MUCH  SHORTER  THAN  THE  NAIVE  NT0Tw*2,  BEING 
GIVEN  BY  THE  FOLLOWING  FORMULA,  DECOMPOSE  MOT  INTO 
2**K2  A 3**K3  A 5A*K5  a ,,,,  LET  SUM2  a 2AK2,  SUMF  a 3aK3  ♦ Saks 
♦ AND  NF  a K3  ♦ K5»  ♦ ,,,,  THE  TIME  TAKEN  BY  A MULTI- 

DIMENSIONAL TRANSFORM  ON  THESE  NTOT  DATA  IS  T « T0  ♦ NTOTa(TIA 

T2ASUM2AT3ASUMrATRANF) , ON  THE  CDC  3300  (FLOATING  POINT  ADD  TiMt 
OF  SIX  MICROSECONDS),  T a 3000  ♦ NTOTa ( B00t<»3ASUM2468ASUMFA 
320ANF)  MICROSECONDS  ON  COMPLEX  DATA,  IN  ADDITION,  THE 
ACCURACY  IS  GREATLY  IMPROVED,  AS  THE  RMS  RELATIVE  ERROR  IS 
BOUNDED  BY  3a2aa (-8) aIUM ( FACTOR ( J) a* 1 , 5) , WHERE  B IS  THE  NUMBER 
OF  BITS  IN  THE  FLOATING  POINT  FRACTION  AND  FACTOR(J)  ARE  THE 
PRIME  FACTORS  OF  NTOT, 

PROGRAM  BY  NORMAN  BRENNER  FROM  THE  BASIC  PROGRAM  BY  CHARLES 
RADER,  RALPH  ALTER  SUGGESTED  THE  IDEA  FOR  THE  DIGIT  REVERSAL, 

MIT  LINCOLN  LABORATORY,  AUGUST  1967,  THIS  IS  THE  FASTEST  AND  MOST 
VERSATILE  VERSION  OF  THE  FFT  KNOWN  TO  THE  AUTHOR,  SHORTER  PRO- 
GRAMS FOURl  and  F0UR2  RESTRICT  DIMENSION  LENGTHS  TO  POWERS  OF  TWO, 
SEE—  IEEE  AUDIO  TRANSACTIONS  (JUNE  1967),  SPECIAL  ISSUE  ON  FFT, 

THE  DISCRETE  FOURIER  TRANSFORM  PLACES  THREE  RESTRICTIONS  UPON  THE 
DATA, 

1,  THE  NUMBER  OF  INPUT  DATA  AND  THE  NUMBER  OF  TRANSFORM  VALUES 
MUST  BE  THE  SAME, 

2.  BOTH  THE  INPUT  DATA  AND  THE  TRANSFORM  VALUES  MUST  REPRESENT 
tOUlSPACEO  POINTS  IN  THEIR  RESPECTIVE  DOMAINS  OF  TIME  AND 
FREQUENCY,  CALLING  THESE  SPACIN68  DELTAT  AND  OELTAF,  IT  MUST  BE 
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C TRUE  THAT  DELT AFs2aP 1/ ( NN ( I ) *DKLT A T ) , OF  COURSE,  DELTAT  NEED  NOT 

C BE  THE  SAME  FOR  EVERY  DIMENSION, 

C 3,  CONCEPTUALLY  AT  LEAST,  THE  INPUT  DATA  AND  THE  TRANSFORM  OUTPUT 

C REPRESENT  SINGLE  CYCLES  OF  PERIODIC  FUNCTIONS, 

C 

C EXAMPLE  I,  THREE-DIMENSIONAL  FORWARD  FOURIER  TRANSFORM  OF  A , 

C COMPLEX  ARRAY  DIMENSIONED  32  BY  25  BY  13  IN  FORTRAN  IV, 

C DIMENSION  OATA( S2,2S, 13},WORK(50),NN(3) 

C COMPLEX  DATA 

C DATA  NN/32,25,13/ 

C DO  I 1*1,32 

C DO  1 J»l,25 

C DO  I K«l,13 

C I DATAd',  J,K)*C0MPLEX  VALUE 
C CALL  FOURT(DATA,NN, 3,-1,  1,1nORK) 

C 

C EXAMPLE  2,  ONE-DIMENSIONAL  FORWARD  TRANSFORM  OF  A REAL  ARRAY  OF 
C length  64  IN  FORTRAN  II, 

C DIMENSION  0ATA(2,6«) 

C DO  2 1*1,64 

C DATACl, I)*REAL  PART 

C 2 DATA(2,I)«0, 

C CALL  F0URT(DATA,64, 1,-1,0,0) 

C 

; DIMENSION  DATAd), NN(1),IFACT(32),N0RK(1) 

C CDC  66B0  INITIALIZATION 

i WR*0, 

1 WI*0, 

> WSTPR*0, 

> NSTPI*0, 

' TWCPI«6, 283185307 

i IF(NOIM-1)920, 1, 1 

) I NT0T*2 

1 DO  2 IDIH«1,N0IM 

I IF(NNdDlM)}920,920,2 

; 2 NTOT*NTOT*NNdDIM) 

C 

C MAIN  LOOP  FOR  EACH  DIMENSION 

C ^ 

i NP1«2 

I DO  910  IDIM*1,N01M 

) NaNNdDIM) 

, NP2*NP1*N 

f IF(N-1)920,900,5 

C 

C FACTOR  N 

C 

J 5 M»N 

) NTWOaNPl 

J IFal 

I IDIV*2 

; 10  IQUOT*M/IOIV 

} ZREMaM-lOIValGUOT 

( ZFdQUOT*IOiy)50,ll,ll 
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: 11 

> 1? 

I 

i 

) 2P 

) iP 

I 

> 

i 31 

» 32 

> 

> 

r 

) 40 

) 

J 50 

51 

I 

1 60 

C 
C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


IK(IRFM}20» 12>2e 

NTnO*NTIn04NTIn0 

M«1QU07 

GO  TO  10 

ID1VS3 

IQU0T»M/1D1V 

IREM»H*IOIV«IOUOT 

IKiaUOT»IOIV)60#3l,31 

IF(1PEM)40,32»40 

UACT(1P)«I0IV 

p«iauo7 

GO  TO  30 
IDIVCI01V42 
GO  TO  30 

ir(IREM)60«51»60 
NTW0«NTK04NTW0 
GO  TO  70 
IKACTCIF)»m 

SEPARATE  FOUR  CASE5»»» 

1,  COMPLEX  TRANSFORM  OR  REAL  TRANSFORM  FOR  THE  4TH,  5TH,ETC, 
DIMENSIONS, 

2,  REAL  TRANSFORM  FOR  THE  2N0  OR  3HD  DIMENSION,  METHOD- 
TRANSFORM  HALF  THE  DATA,  SUPPLYING  THE  OTHER  HALF  BY  CON- 
JUGATE SYMMETRY, 

3,  REAL  TRANSFORM  FOR  THE  1ST  DIMENSION,  N ODD,  METHOD-* 
TRANSFORM  HALF  THE  DATA  AT  EACH  STAGE,  SUPPLYING  THE  OTHER 

half  by  conjugate  symmetry, 

4,  REAL  transform  FOR  THE  IST  DIMENSION,  N EVEN,  METHOD-* 
TRANSFORM  A COMPLEX  ARRAY  OF  LENGTH  N/2  WHOSE  REAL  PARTS 


C ARE  THE  EVEN  NUMBERED  REAL  VALUES  AND  WHOSE  IMAGINARY  PARTS 

C ARE  THE  ODD  NUMBERED  REAL  VALUES,  SEPARATE  AND  SUPPLY 

C THE  SECOND  HALF  BY  CONJUGATE  SYMMETR'', 

C 

I 70  N0N2«NPiw(NP2/NTH0) 

> ICASE*! 

» IF(IOIM.4)71»90,90 

' 71  IFCIFORM)72,72,90 

) 72  lCAie»2 

) IF(IDlM.n73#73,90 

1 73  1CA8E«3 

JF(NTWO*NPn90,90,74 

> 74  ICASE«4 

( NTW0«NTW0/2 

I NaN/2 

i NP2aNP2/2 

I NTOTbNTOT/2 

' Ia3 

1 DO  60  Ja2,NT0T 

) OATA(J}aDATA(I) 

) 80  lalt2 

90  IlRNGaNPl 

! IF(ieA8E*2)100,95,100 

i 95  IlRN6aNP0*(i4NPREV/2) 


I 
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C 

C SHUFFLE  ON  THE  FACTORS  OF  TWO  IN  N,  AS  THE  SHUFFLING 

C CAN  BE  DONE  BY  SIMPLE  INTERCHANGE,  NO  WORKING  ARRAY  IS  NEEDED 

C 

t too  IF (NTWO*NP1)6P0,&00, IIB 

i 110  NR2HF»NP2/? 

jBl 

DO  150  I2«1,NP2,N0N2 
IF(J«I2) 120,130,130 
IlMAXaI24'N0N2«>2 
DO  125  I1bI2,I1MAX,2 
DO  125  I3b11,NT0T,NP2 
J3aJ«I3*I2 
TEmPRbDATA(I3) 

TEMPIaOATA(X3^1) 

DATA(I3)aOATA(J3) 

DATA(I3tl)aDATA(J3«n 
DATA(J3)aTEMPR 
DATA(J34naTEMPI 
HaNP2HF 

IF(J»M)150,150,lit5 
JaJ«H 
MaM/2 

IF(M.N0N2)  150,  H0,  l«0 
JaJ^M 

MAIN  LOOP  FOR  FACTORS  OF  TWO,  PERFORM  FOURIER  TRANSFORMS  OF 
LENGTH  FOUR,  WITH  ONE  OF  LENGTH  TWO  IF  NEEDED,  THE  TWIDDLE  FACTOR 
WaEXPCISlCNB2*PI*S0RT(-l)*M/(4*MMAX)),  CHECK  FOR  WalS IGN*SORT ( • 1 ) 
AND  REPEAT  FOR  WalS IGN*SQRT (• 1 ) wCON JUGATE ( W ) , 

N0N2TaN0N2^N0N2 
IPARaNTWO/NPl 
1F(IPAR*2) 350,330,320 
IPARaIPAR/4 
GO  TO  310 

DO  340  Il8l,IlRNG,2 
DO  340  J3all,N0N2,NPl 

DO  340  KlaJ3,NT0T,‘0N2T 

• K2aKl4N0N2 

( TEMPRaDATACK2) 

; TEMPIaOATA{K2^1) 

• DATA(K2)aDATA(Kl)«TEMPR 

' DATA(K24l)a0ATA(KUl).TEMPI 

t DATA(Kl)aDATA(Kl)tTEMPR 

> 340  OATA(KUl)aOATA(Kl4l)4TEMPI 

) 350  HHAXaN0N2 

360  IF(MMAX«NP2HP)370,6e0,600 

! 370  LMAXaMAX0(NON2T,MMAX/2) 

i IF(MMAX«NON2}405,405,3S0 

I 380  THETAa«TW0PI*FL0AT(N0N2)/FL0AT(4*MHAX) 

i 1F(I8ICN)«00,390,390 

, 390  THETAa.TMETA 

' 400  WRaCOSITHETA) 


t 

\ 

t 120 
1 


> 


» 125 

t 130 

1 140 

145 

» 

> 

I 150 

C 
C 
C 
C 
C 

c 


' 310 

i 320 

) 

) 330 

k 
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> WlaSINCTHETA) 

) MSTPRao2,*K^I*wl 

1 WSTPI«2,*NR*WI 

. 4P5  DO  570  L»N0N2f LMAX»N0N?T 

I m«l 

i IF(MMAA«NON2)420,420»410 

I 410  M2RbMR*WR«WI*mI 

> W2la2,*WR*wi 

> w3Rai«i2R*WR«M21*Ml 

' N3lBW2R*MI«W2I*WR 

) 420  DO  530  Ilal>IlRN6,? 

> DO  530  J3«I1»N0N2,NP1 

5 KHlNaJ3tIPAR*P 

I IF (MMAX-NON2)430#430»44O 

> 430  KMINaJ3 

5 440  KDIFaIPAR*MMAX 

\ 450  K8TEPa4*KDIF 

i DO  520  KlaKMlN,NITOT,KSTEP 

> K2aKl4KDIF 

f K3aK2^KDIF 

J K4aK3+KDIF 

> IF(HMAX.NON2}460#460>4e0 

) 460  UlRaDATA(Kl}^DATA(K2) 

UlIaDATA(Kl+l)40ATA(K24l) 

» U2R80ATA(K3)^DATA(K4) 

i U2laOATA(K3fl)+DATA(K4+l) 

I U3RaOATA<Kl)*DATA(K2) 

i U3IaDATA(Kl+l)*DATA(K24l) 

. IF(ISI6n)470»475#475 

' 470  U4RaDATA(K3fl)-DATA(K4tl) 

J U4laDATA(K4)»DATA(K3) 

> GO  TO  510 

) 475  U4RaDATA(K4tl}*0ATA(K3f 1) 

U4laOATA(K3)«DATA(K4) 

> GO  TO  510 

i 480  T2RaW2R*DATA(K2)*W2laDATA(K2>l) 

I T2laW2R*0ATA(K24l)*M2I*DATA(K2} 

> T3RaWR*DATA(K3)«Wl*DATA(K3tl} 

. T3Ial«(R*DATA(K3fl)4WI*DATA(K3) 

' T4Raw3RaDATA(K4)-W3I*DATA(K4*l) 

1 T4laW3R*DATACK4^1)*W3I*DATA(K4) 

> UlRaOATA(Xl)+T2R 

) UllaDATA(Kl4l)fT2I 

i U2RaT3R4T4R 

» U2laT31*T4l 

I U3RaDATA(Kl)«T2R 

\ U3laDATA(Klfl)«T2I 

> IP(X8ION)490, 500^500 

> 490  U4RaT3Z»T4I 

' U4laT4R«T3R 

\ 60  TO  510 

f ( 500  U4RaT4l«T3I 
) U4laT3R«T4R 

510  DATA(Kl)auiR4U2R 
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DATA(Ki^n»UlI  + U2I 
i DATA(K2)«U1R+U4R 

> DATA(K2tnaU3l4U4I 

i DATA{K3)»U1R<»U2R 

) DATA(K3*l)»UH-U2I 

' DATA(K4)«U3R«U4R 

< 520  DATACK4^1)aU3I*U4I 

> KMIN«4*(KMIN«J3)+J3 

> KDlFaKSTEP 
IKKOIF-NP2)450,530,530 

! 530  CONTINUE 

i MaMMAXaM 

I IF(I8IGN)540,550,550 

> 540  TEMPRaWR 

> WRaaNI 

’ WlaaXEMPR 

) GO  TO  560 

> 550  TEMPRaWR 

) WRbnI 

WlaTEMPR 

! 560  IK(M«LPAX)565,565,410 

i 565  TEMPRbwR 

I WRakxRaWSTPRaWIaWSTPI  + WR 

i 570  WlaWlaNSTPRfTEMPRawSTPI+WI 

. IPARa3aIPAR 

' nnaxbmmaxxmmax 

» GO  TO  360 

C 

C MAIN  UOOP  FOR  FACTORS  NOT  fcQUAU  TO  TWO,  APPLY  THF  TWIDDLE  FACTOR 
C W«EXP(I8IGN*2aPIwS0RT(al)*(J2-l)*(JlaJ2)/(NP2*IFPl))»  THEN 
C PFRFORM  A FOURIER  TRANSFORM  OF  LENGTH  IFACTC1F)»  MAKING  USE  OF 

C CONJUGATE  SYMMETRIES, 

C 

) 600  IF(NTWOaNP2}605f 700^700 

I 605  IFP1«N0N2  ^ 

iFal 

! NPlHFaNPl/2 

I 610  IFP2aIFPl/IFACT(IF) 

I JlRNGaNPa 

> IF(lCA8Ea3)612f611f 612 

> 611  JlRNGa(NP2*IFPl)/2 

' J28TPaNP2/irACT(IF) 

« JlRG2#(J2S7PxIFP2)/2 

> 612  J2MlNaUIFP2 

) IF(IFPlaNP2)615r640,640 

; 615  DO  635  J2aJ2MIN,IFPl,IFP2 

> THETAaaTW0PIaFL0AT(J2al)/FL0ATCNP2) 

i IF(ISIGN}625#620,620 

I 620  THETAaaTHETA 

; 625  8INTHa8IN(THETA/2,) 

> WSTPR«a2,a8INTH*8INTH 

' W8TPI«81NCTHETA) 

\ WRaWSTPRtl, 

) WlaWSTPI 
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' JlMIN«J2tUPl 

DO  635  JiRNG, IFPl 

’ I1MAX»JHI1RNG*2 

> 00  630 

( 00  630  I3»IlrNT0T»NP2 

> J3MAXb134IPP2«NP1 

> 00  630  J3bI3» J3MAX,NP1 

' TE.MPR«DATA{  J3) 

I 0ATA(J3)»0ATA(J3)*WR«0ATA(J34i)BWI 

> 63U  DATA(  J3^n»TtMPR*lill<.0ATA(J3+n*WR 

) TEMPRbNR 

MRbNRAWSTPRbHIamSTPUmR 
• 635  WIsTFMPR*i«3TPI^WI*W8TPR^WI 

i 6R0  THETAb-Th^OPI/PLCAKIKACTCIK)) 

I IF(I8IGN)650,6R5,645 

» 6<I5  THETAb-THETA 

. 650  SINTHbSIN(THETA/2,) 

' WSTPR»*2,*3INTH*SINTH 

I WSTPI»8IN(THfcTA) 

> K8TFPk2*N/IFACT(IF) 

» KRANC»KSTtPB(IFACT(IF)/2)+l 

00  696  Ilsl,URNC,2 
; 00  698  I3«I1»NT0T,NP2 

i 00  690  KMIN«1,KHANG,KSTEP 

I JlMAX«I3+JiRN0»IFPi 

> 00  680  JlsI3f JIHAX, IFPl 

. J3MAX*JI*IFP2»NP1 

’ 00  680  J38J1, J3MAX,NP1 

J J2MAX«J3+IFPi«.lFP2 

I K«KMIN4(J3-J14(J1-I3)/IFACT(IF))/NPIHF 

) IF(KMIN«n655,655f  665 

655  SUMR«0, 

! SUHIB0, 

} 00  660  J2bJ3, J2MAX, IFP2 

I SUMRBSUMRtDATA(J2) 

» 660  SUMIbSUMUOATAC  J2*n 

> W0RK(K)b8UMR 

' W0RK(K*1)b8UMI 

I GO  TO  680 

f 665  KC0NjBKt2B(N«iKMIN4-l) 

) J2BJ2MAX 

8UMRbDATA(J2) 

; SUMIbDATA(J241) 

i Ol,OSRB0, 

I OLO81B0, 

5 J2BJ2-IFP2 

. 670  TEMPRB8UMR 

' TEMPIbSUMI 

) 8UPRbTM0MR*8UMR«0L0SR«0ATA( J2) 

> SURlBTM0WR*SUMI«0L0SIt0ATA(J2tn 

) 0LD8RBTEMPR 

OLDSIbTEMPI 
! J2BJ2*1FP2 

) IF(J2-J3)675,675,670 

I 675  TEMPRbWR*8UMR.0LDSR40ATA(J2) 
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7EHPI«WI*SUMl 

I»0RK(K)»TEMPR-TEMP1 

NORK(KCONJ)aTEMPR+TCMPI 

TEMPRsWRaSUHI-OLDSI^DATAC  J2tn 
TFPP1»NI»SUMR 
WORK(K41)*rtMPR4TEMPI 
NORK(KCONJfn»TEHPR*TEMPI 

6eu  CONTINUE 

IKKHIN«n68‘j,685,6e6 
68b  WRanJSTPRfl , 

WlrWSTPI 
GO  TO  690 
686  TfcMPRaWR 

(»R»WR*N8TPR-WI**«iSTPIfNR 
NlaTEMPRawSTPUWIaWSTPR  + WI 

690  T*<0WR«WR4WR 
ir(ICASF-3)69?,69l,692 

691  IKIFPl-NP2)69b,692,692 

692  K»l 
I2MAXsI3*NP2«NPl 

00  693  I2«I3, I2MAX,NP1 
DATA(I2)xwORK(K) 

DATA(I2fl)«H0RK(Ktl) 

693  KaKA2 

GO  TO  698 
C 

C COMPLETE  A REAL  TRANSFORM  IN  THE  1ST  DIMENSION,  N 000,  BY  CON- 

C JUGATE  SYMMETRIFS  AT  EACH  STAGE, 

C 

695  J3MAXaI3fIFP2-NPl 

DO  697  J38I3, J3MAX,NP1 

J2MAX«J3+NP2*J28TP 

DO  697  J2sJ3, J2MAX, J28TP 

JlMAXaJ2fJlRC2»IFP2 

J1CNJ«J3+J2MAX^J2STP»J2 

DO  697  JlaJ2, JIMAX, IFP2 

K«1+J1*I3 

DATA(  JDaWORKCK) 

DATA(  JUn«WORK(K^l) 
ir(Jl*J2}697r697,696 

696  DATA( J1CNJ)«W0RK(K) 

DATA (JlCNJ+1) ■•WORK (Ktl) 

697  JICNJ«JICNJ«IFP2 

698  CONTINUE 
IFbIF+1 
IFPl«IFP2 

IF(IPPI-NP1)700»700»610 

C 

C COMPLETE  A REAL  TRANSFORM  IN  THE  1ST  DIMENSION,  N EVEN,  BY  CON. 
C JUGATE  SYMMETRIES, 

C 

700  60  TO  (900,800*900,701), ICA8F 

701  NHALFbN 
NbNTN 
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THETA««TW0PI/FL0AT(N) 

’ IF(IS1GN)703»  70i?,  702 

: 702  THETAa»THtTA 

1 703  SINTH«8IN(THtTA/2,) 

. WSTPR«*2,*SINTH*SINTH 

. iNSTPlaSIN(THtTA) 

’ WR«WSTPR+1, 

i WlaNSTPI 

I IMiNs3 

) JMIN«2*NHAL7»1 

GO  TO  725 
• 710  JaJMIN 

i DO  720  l«IMIN,NTOTrNP2 

» SUMR«CDATA(I)+DA1 A(J))/2, 

> SUMIa(DATA(l-H)+DATA(J^n  )/2, 

. DIFRs(DATACI)-»DATA(J))/2, 

’ DIH»(DATA(Ifn*DATA(  J*n  )/2, 

t Ttf<PRsWR*SOHUWI*DIFR 

> tfmpi«wiasumi*wr*difr 

» DATA(I)«SUMR*TEMPR 

DATA(I+n*DIFI  + TEMPI 

: data(J)*sumr*tempr 

J DATA( J^lJs.DIFI+TEMPI 

I 720  J»JtNP2 

; IMINaIMIN+2 

. JMlNsJMIN-2 

' TEMPRaWR 

1 WRbwR*WSTPR*WI*'^STPI4-WR 

I >ilBTEMPRot^3TPI  + wI*»<STPR  + wi 

) 725  IF(IMIN-JM1N)710»730,7A0 

730  IF(IS1GN)731,7R0»7«0 

; 731  DO  735  IbIMIN,NT0T,NP2 

735  DATA(1A1)b-DATA(UI) 

740  NP2BNP24NP2 

NTOTBNTOTtNTOT 
JbnTOT+1 
' IMAXBNT0T/2tl 

I 745  IMINbIHAX»2*NHALF 

} IBIHIN 

> CO  TO  755 

750  DATA(J)bDATA(I) 

: DATA(Jfl)B.OATA(Ul) 

» 755  lBlf2 

I J«J-2 

) IF(I*IMAX)750»760,760 

. 760  DATA(J)BOATA(IMIN)-OA7A<lMINfl) 


DATA(J+1)b0, 

IF(I«J)770»780»780 
765  DATA( J)bDATA(I) 

DATA(J+l)BDATA(I^n 
770  I«I»2 

J«J»2 

IF(Ib1MIN)775,775,765 
775  DATA(J)BOATA(IMlN)*DATA(IPINtn 
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DATA(  J+ns0, 

IMAX»IMIN 
GO  TO  7/15 

7B0  DATA(1)bDATA(1)+DATA(2) 

DATA(2)s0, 

GO  TO  900 
C 

C COMPLETE  A REAL  TRANSFORM  FOR  THE  2ND  OR  iRO  DIMENSION  BY 

C CONJUGATE  SYMMETRIES, 

C 

800  IF(1  lRNG*NPn805,900,900 
805  DO  860  I3sl,NT0T,NP2 
I2MAXaI3+NP2»NPl 
DO  860  12*13, I2MAX,NPl 
IMINsI2^IlRNG 
IMAX«I2^NP1«2 
JMAXs?*I3+NP1.IMIN 
IF( 12*13)820,820,610 
810  JMAXaJMAX+NP2 
820  IFdDIM. 2)850, 850, 830 

830  JaJMAX+NP0 

DO  0a0  lalMIN, IMAX,2 
OATA(I)*DATA(J) 

DATA( I41)b.DATA( J+l) 

890  JaJ*2 
850  JaJMAX 

DO  860  I»IMIN, IMAX,NP0 
DATA(l)aDATA(J) 

DATA(I+l)a*DATA(j4l) 

860  JaJ*NP0 
C 

C END  OF  LOOP  ON  EACH  DIMENSION 

C 

900  NP0*NP1 
NPlaNP2 
910  NPREVaN 
920  RETURN 
END 

ROUTINES  CALLEDI 

MAX0  , FLOAT  , COS  , SIN 

OPTIONS  a/0N,/0P|3 

BLOCK  LENGTH 

FOURT  3972  (015990)* 

**C0MPILER  CORE** 

PHASE  USED  free: 

DECLARATIVES  00622  19690 
executables  01583  13729 
ASSEMBLY  02935  17517 
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